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Part I 


PROBLEM DEFINITION AND 
CLASSICAL SOLUTIONS 


Chapter 1 


STIFF? WHAT IS IT? 


1.1 The initial value problem 


A stiff problem is a particular case of an initial value problem. We define the CAUCHY’s 
problem or initial value problem a system with an arbitrary number of ordinary dif- 
ferential equations of first order endowed with a given starting value: 


m = — axtxb, (1.1) 


In this case, there are N equations and y : [a,b] — RY, f : [a,b] x RN > RN, 
Ya € RN. 

The step by step numerical algorithm builds an approximate solution by an inter- 
ative procedure starting from t = a. Note y(t,) the exact solution at the point t» of 
the associated interval division 


a = İQ <: <tn <- CİM =), 


and yn, the approximate value computed with the numerical method. 
For the unicity of the exact solution we put the following conditions: 


e the system function f is continous on 
D = {(t,y) |t € [a,b], İyi) — yil < d}, 
where d € R4; 
ə f satisfies the LIPSCHITZ’s inequality 
MoNS Luv], arx, WED (Lə) 
For almost all differential equation systems, the value L(b— a) is of hundreds order. 
For such systems the classical step by step methods, like RUNGE-KUTTA’s processes 
or ADAMS’s formulae, are satisfactorily (in the sense of a small error). 


Unfortunately, there are some exceptions, for instance, the case when the function 
variation is very strong. When the value L(b — a) exceeds the thousands order, a 
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variety of restrictions are imposed to the classical methods, especially on the stepsize, 
so there are useless. Such a system we classified to be of stiff kind one. 


Remark. Applying (1.2), we deduce the existence of a value A > 0 so that the 
inequality 
lut + A) — vt + A)İ s lut) = vll, vt e [a,b], (1.3) 
holds for two arbitrary exact solutions of the differential system, u and v, which are 
starting from different initial values. The value A depends on L. If the value L is 
large, the value A is very small. The inequality (1.3) is important for the numerical 
integration because it tell us that the approximate solution can follows an other exact 
solution which start from an initial value which is closed to that of the exact solution 
and the approximation error (the difference between the approximate solution and the 
exact solution) can be controlled. 


A more general condition is that of the dissipativity of the system function: 
< f(t,u) — ft, v), u — v >< ullu — oll”, (1.4) 


where u is the constant of dissipativity. When y, f(t, y) € CN, the dissipativity con- 
dition can be reformulated: 


Re < f(t,u) — f(t, v), u — v >< İllu — oll”. (1.5) 


A LIPSCHITZ continous function satisfies the dissipativity condition, but the reverse 
implication is not valid. The value u can be lowest than the value L, also a negative 
one. 


Remark. The dissipativity condition confirms the existence of a value A for which 
lult + A) — v(t + Al e lut ll, Ve € fa, b]. (1.6) 


For u <0, the difference between two arbitrary solution is a decreasing function. 


1.2 The stability of the exact solutions 


Let y’ = Ay, A € R. The solutions family of this equation, can be represented by 
{y| y(t) = ceX, c€ R). Note that, when \ > 0, for some large t values, the difference 
between two arbitrary solutions increases exponentially, even if the difference between 
the initial values is very small. In such a case, applying an approximate method, 
the results are unforeseeables. If A < 0, the difference between two arbitrary exact 
solutions decreases when t — oo. These remarks can be extended for the case 


y =A, yEC (1.7) 


(discussion on Re à > 0 or ReA < 0). An equation of (1.7) form is referred to as the 
scalar test equation. 

The scalar test equation is used in the numerical stability study of the step by step 
methods, especially for the multistep formulae. 

Let f : [to, co) x RY > R”, y(-, to, go) the exact solution of the CAUCHY’s problem 
with y(to) = yo, and y an arbitrary solution for the differential system. 
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Definition 1.1. The solution g is stable in LIAPUNOV’s sense if, for any € > 0, there 
is a value d(¢,to) such that 


lC) — y(t, tə, yo)l| ce, VE € [to, 00), (1.8) 


holds for any yo satisfying ||yo — y(to)|| < ö(e, fo). The solution g is asympotically 
stable if, for any € > 0, there is a value ö(fo) such that (1.8) holds for any yo satisfying 
Ivo — e(to)İİ < (to). 

In the particular case of a linear system, y’ = Ay, if the solution y = 0 is stable 
(asymptoticallly stable), all the system solutions are stable (asymptotically stable). 
The solution y = 0 is asymptotically stable if all the eigenvalues of the matrix A have 
some negative real parts, that means that the characteristic polynomial of A is an 
HURWITZ polynomial. 

The stability of the exact solutions is out of the present study subject. 

The following supposition is esentially in our study. 


The exact solutions of the numerical solved problems with initial 
values are asymptotically stable 


By this supposition we eliminate the possibility of a numerical instability due to the 
instability of the exact solutions. 


1.3 The stability of the numerical solutions 


If the error propagation mechanism is under control, the numerical method is stable. 
In mathematical terms, this fact can be expressed in the following form. 


Definition 1.2. A method (formula, process, scheme) is stable if, for each differential 
equation with all solutions asymptotically stable ones, there are two values, K and ho, 
so that 

lyn — nll < Kilo — Voll, Vh: O<h<ho, (1.9) 


holds, where y, and 7, are the values of the approximate solutions in the nüh point 
of the interval division, computed with the stepsize h and starting from the value 
Yo, respectively Yo. The method is asymptotically stable if, for a given stepsize, the 
perturbations in the numerical solution do not grow from a step to another, i.e. 


lyn = Yall < İli = Gr-al- (1.10) 


Remark. The numerical stability in the sense of definition 1.2, tell us that a small 
perturbation produced to one step of integration determines only a small perturbation 
in the numerical solution for all future points, when we use a small stepsize. 
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1.4 Order and consistency 


Applying the numerical method, we find out some approximate values y, of the exact 
values y(t,). Ideally, yn = y(tn), ie. the numerical algorithm is exact. Unfortunately, 
this request can be satisfied only in the case when we use a multistep formula to 
integrate an equation with a polynomial solution. Generally, the difference between 
the approximate solution and the exact solution must vanish when the division points 
number becomes infinitely. That means that the approximation can be improved 
growing the number of the division points, i.e. increasing the computational effort. The 
consistency express this request in the theoretical case when the numerical algorithm 
can be applied without truncation errors. 

The method order is the convergence order of the approximate values set at h > 0 

Example. Let the onestep method 


Ynt1 = Yn + hÖ(t,, Yn, h). (1.11) 
Note 
Zn+1 :— y(t) + haz? (tn, y(t), hazı) : TE = y(tn41) — Zn-. 


The value TE is the local discretization error. The global error can be expressed by 
the difference 


Y(tn41) — azi — TE + (2n41 — yası) 


If 
TE = ChPtlo(t, y(t,)) + O (h”12), 


where y depends on the system function (for almost all multistep formulae, is a p 
order derivative of the system function), then the algebraic condition of consistency is 
p > 1. The p value is the accuracy order of the numerical method. 

The notion of accuraccy supposes the consistency condition and a small error con- 
stant C, smallest than the stepsize A. 


1.5 Behaviour of classical methods for the stiff 
case 


We study a suggestive example. Let 
yz Ay, te[0,7], (0) = yo 


and a division of the integration interval with the constant stepsize h > 0 and the 
points tn = nh, n > 0. The exact solution y satisfies the reccurence relation y(tn41) = 
e“"y(t,,). The most simple numerical procedure, for finding an approximate value y, of 
the exact value y(t,), is the EULER’s explicit rule. For the given system, the iterative 
process is 

naq. Yn + hAyn, n>0 
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Set S(z) = e, R(z) = 1 +z. Then the reccursive equations of the solutions are the 
followings: 


y(tn41) = S(hA)y(tn), Yn+1 = R(hA)yn- 


The global error at the point tn is the value ep = y(tn)—yn. The reccursive equation 
is en41 = R(hA)e,+T(hA)yn, where T = S—R is the truncation function. The solution 
of the reccursive equation starting from e9 = 0 is en4, = io RI(hA)T(hA)Yn—j. 
Thus, the global error can be approximated by 


J . 
lenll Sm max RUAJ max İ70:4)yil. 


If the numerical solution is stable, i.e. || R(RA)|| < 1, and is of p order of consistency, i.e. 
| (hA)y|| = O(h?**)||y||, then an approximation of the global error is İle, li = O(R?) 
and the method converges, i.e. e, — 0, when n — oo. For the EULER’s explicit rule 
p= İ, since 
Hö... 
where A, are the eigenvalues of the system matrix, A. Thus we can explain the use of 
the consistency order as a measurement of the convergence order of the method. 
The stability condition || R(hA)]| < 1 means that 


İH EBA, 71.....N. 


The set 
A= {z||R(z)|| x 1. 


is the absolute stability domain. For the EULER’s explicit rule, R(z) = 1 + z, and 
the stability domain is the set of complex values which are inside to the unitary circle 
{z € Cİİl 4 zl = 1}. 

For usual differential equations, |Amaz| = maxi<i<n |A;| is not a very large value 
and the restriction imposed by the stability condition, z; = hà; € A, i = 1...., N, is 
satisfied for some reasonable stepsizes. A classical example of stiff system is the one 
for which |Amax| is very large. In this case, neither the stability condition, nor the 
error boundness can not be satisfied without the strong restriction of the stepsize, in 
fact introducing a large computational overhead. The term ”large” is relative to the 
value 1/7. Indeed, also an equation with a |Amax| small value can be classified as a 
stiff one if the requested solution is defined on a very large integration interval. 


Remark. Let the CAUCHY”s problem associated to the nonlinear equation y’ = 
f(t,y) with o, a numerical solution. In a y(t) neighbourhood we can approximate the 
derivative 


y(t) s F(t, (4) + IEU) - et). 
where J(t) is the JACOBIAN’s matrix (Of(t,y)/Oy)(t, yp(t)). Then the global error 
e(t) = y(t) — y(t) can be approximated by 
e(t) x J(t)e(t). 


Thus, the error evolution depends on the eigenvalues of the matrix J (t). 


14 CHAPTER 1. STIFF? WHAT IS TT? 


A numerical method for solving stiff problems is a method for which the stability 
condition does not ask that |hA;| remains small also when the corresponding com- 
ponent in the exact solution is negligible. Practically, there are such methods. For 
example, the EULER’s implicit rule, with the reccursive relation 


Untl = Yn + h faa, (1.12) 


where fayı = f(tnt1,Yn41), has the requested property. The stability domain is the 
set of the complex points lying outside the unitary circle {z € CİİL— z| = 1}. 

The stepsize restriction is only one characteristic of the stiff phenomenon. 

The numerical methods based on polynomial interpolation are succesfully in the 
integration of a stiff system only if they are implicit, and, consequently, suppose some 
procedures for solving the nonlinear difference equations. Solving implicit nonlinear 
equations associated to stiff problem, the use of the method of simple iterations is 
out of question. For example, for the formula (1.12), the method of simple iterations 
concerns the iterative process 


məebl m 
or = Un thf (tori yl) . m? 
where gə, is given by an explicit method. The convergence condition of this iterativ 
process 1s 


o 
hp (Hitnu) < l, 


which can not be satisfied for a reasonable value 2) relative to the integration interval 
length. The most popular method for solving the implicit equations, in the stiff case, 
is the NEVVTON”s iterative process. "The JACOBIAN’s matrix must be computed to 
each iterative step, which is very expensive in the sense of the computational effort. 


1.6 Definition 


The essence of the stiff phenomenon consists in the fact that the exact solution includes 
some components with a very fast decreasement that can be very hard to be followed 
by the numerical solution given by a step by step iterative process. 

We mention some pragmatic definitions: 


ə The stiff equations are the equations for which some implicit methods work 
better than the explicit ones. 


e A problem is stiff in a given integration interval if, for a given numerical code, 
the stepsize must be very strongly reduced. 


e The stiff differential equations are wrong conditioned in the computational sense. 


e An ordinary differential equation system is stiff in a given integration interval 
(0, T] if, in the exact solution, there is at least one component with a very large 
variation relative to the value 1/7. 


1.6. DEFINITION 15 


Example 1. Let the scalar equation 
y(t) = f(t,y) = Aay HEO AF, 120, A€ 0. 
The exact solution of the associated CAUCHY’s problem is 
y(t) = F(t) +e“ [yo — F(0)], 


We notice that, after a short time, the transitory component e™[yo— F(0)), also referred 
to as the stiff component, has no more a signifiant influence on the solution value. 
There is a subinterval in which the slow varying component, F(t), also referred to 
as the untransitory component or the smooth component, is prevalent in the exact 
solution. This subinterval covers almost all the integration interval. Applying the 
EULER’s explicit rule with stepsize A to this problem, we get 


Ynti = Yn + hf (tas Yn) = (1: həl — (tn) + En) + AE Ua). 


If the numerical model is correct, the difference yn — F(t,) must decrease: this is true 
iff —2 < hà < 0. Thus, the stepsize is retricted by some numerical stability condition. 

In the transitory stage a problem can not be classified to be stiff one, only when 
the stiff component is negligible. An equation can not be referred to as stiff one (only 
the CAUCHY’s problem associated to it), since the term is relative to the initial value, 
the integration interval and the error tolerance. 


Example 2. Some components with different dials of variation can be meet espe- 
cially in the case of the linear systems of the following form 


iti) = Ayi) + FU). (1.13) 


where A is a constant matrix with different dials values of the eigenvalues A,. The 
problem (1.13) is stiff if (FATUNLA[43]) 


i) J: Re(A,) € 0, 

il) 49: Aj] < JA; 

ili) Ak: Re(Az) > 0; 

iv) Am: |Im(\,,)| >> 0, only if Re (A) € 0. 
We define the stiff ratio of the system (1.13) by 


( 
( 
( 
( 


S(A) = max [Re(Aj)|/ min |Re(Aj)]- 
The system (1.13) is stiff if Re (A;) < 0, ? =1,...,N, and S(A) > 1 (CASHİ25)). 


Example 3. An important class of stiff problems with very large dimensions is 
derived from some partial differential equations. For example, let the problem posed 


in [53] 


Ut = Ure, E> 0, OS a Xİ, u(t,0) = u(t, 1) =0, t>0, u(Ü, z) = y(x) 
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and u(t) the approximate values of the exact values u(t, x;), produced by replacing 


Une With 
b7u(t, x; ult, ai +h) — 2u(t,7;) + u(t,a2;-—h 
“ )_ ul ) a ) tul ) 5 
where h is the stepsize on x axis, h = 1/m, m € N. Then u(t), i = 1,...,m — 1 are 


the components of the exact solution of the CAUCHY s linear problem 


—2 1 


y(t) — Ay(t), t>0, A=m? 


The eigenvalues A, of the matrix A, 


752571 1—1.....m —l, 
m 


are uniformly distributed in the interval (—4m5, 0). When m increases (the division is 
more dense) the stiff ratio also increases. 


Example 4. The singular perturbation problems are particular cases of stiff pro- 
blems. The generic form of a such problem is [53] 


zı ə 
.— "m ? 
eyi = g(t, x,y, 6,) g(t.x,y,0) = 0 


In fact this system is a stiff one. For example, let f = g = x+y. Then the eigenvalues 
of the system are e7! + O(1) and —1 + O(c). When e decreases, the stiff ratio is 
increasing. 

Such systems can be solved by some pseudo-stationary methods. We consider the 
reduced problem 


x! = f(t, £,y, 0), 
0 = g(t, z, y,0), 
x(0) =e. 


Unfortunately, there are some systems which can not be simplified in this manner since 
the presence of the control value e determines the evolution of the exact solution. 


Example 5. If the system presents a variation on t, like 


the affiliation to the stiff class is not easy to be establish. The problem is only local 
a stiff one. The variation of the eigenvalues of the JACOBIAN’s matrix can not inform 
us about the perturbation propagation in the exact solution. In practice, a system is 
classified to be stiff one if, at least at a moment (at one £ from the integration interval), 
the associated initial value problem is a locally stiff problem. 
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Table 1.1: Subclasses of stiff problems for (b — a)llg el = 1 


Stiff at limit S = O(10) 
Midle stiffly S = O(107) 


Strong stiffly O(10°) x S < O(10°) 
Extreme stiffly 
Pathological stiffly 


In physics and chemestry, any real system which is modelated by an ordinary 
differential equation system with components supposing some evolutions in different 
dials of time, also referred to as time constants, is treated as a stiff system. A large 
time constant determines the stiff character. The components of practical interest are 
those with a slow variation. Integrating the problem with an explicit method, the 
stepsize must be of small constant order. 


Example 6. A nonlinear system is stiff if at least for one point t, € İla, öl the 
associated linear system y’(t) = (Of /Oy)(ts, y(ts) y(t), t > t, is stiff. 

The meaning of the sign >> in the above mentionated inequalities is relative to the 
integration interval length and the error tolerance. An initial value problem can not 
be classified as stiff one, also when the eigenvalues are of very different values, if the 
integration interval can be compared with the transitory stage, or if the error tolerance 
€ is not very small. 


Definition 1.3. The stiff ratio of the system (1.13) is 
1 
S(A,a,b,£)= max, Al — alg =, (1.14) 


where a, b are the extreme values of the integration interval and e is the error tolerance. 


Depending on the stiff ratio the problems can be classified like in Table 1.1. 
In chemestry dynamics we can often meet a stiff ratio of order O(10°). 


Definition 1.4. An ordinary differential equation system y’ = f(t,y), numerical 
integrated on the interval [a,b], starting from the initial value y(a) = ya, is stiff if the 
function f has a LIPSCHITZ’s constant very large relative to the interval length, and 
the error tolerance. 


Remark. For a system which satisfies only the condition of dissipativity with the 
constant u, the stiff or nonstiff character is decided comparing the value u, the interval 
length, and the error tolerance. Note that the dissipativity constant satisfies the 
inequality 


u > sup w[(Of/Ot)(t, y)] (1.15) 


yeRN 
where glAİ = suppers € z, Az > / € 2,2 is the logarithmic norm of the matrix A. 
The principal characteristics of the stiff problems are ( ATTKENİLİ): 
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e the exact solutions are stable in the sense that small perturbations in the initial 
values are followed only by small perturbations in the exact solutions; 


e trying to solve the problem with the standard methods we get some strict re- 
strictions on the stepsize from the stability conditions. 


A problem necessarly must have these two characteristics to be declared a stiff one. 


Chapter 2 


LINEAR MULTISTEP 
FORMULAE 


The linear multistep formulae are the first methods which were proposed for solving 
the stiff problems (Curtis, 1952). Due to the stiff character, this class of formula 
suffers a variety of restrictions. This fact explains the development of some new class 
of methods, like the nonlinear multistep methods, which are related to the multistep 
linear formulae. 

The fundamentals of the modern theory of stability have been establish by 
Dahlquist (1956). From 1971, when Gear has published his test results on some special 
linear multistep methods, these are often used in the stiff computation. 


2.1 Formula 


Let (to, t1,..- ta) a division of the integration interval with constant stepsize h = 
tig, — ti 1 = 0,..., M — 1. A linear multistep formula is represented by an iterative 
process of the followings form: 


k k 
ə Ak—iYn-i — h ə Brifn—i = 9, 
0 0 


(2.1) 


where fan-i = f(tn—i,Yn-i) and M > n > k. If (ağ + 92)(az + 92) 2£ 0, then k is the 
formula stepnumber. 

The formula is an explicit one if 3, = 0, else it is an implicit one. If a, = 0, Bk 4 0, 
then the formula is an extended one. 

At the step n, knowing the (k — 1)'" previous approximate values, we get the kt 
value by writing the formula in the form 


AkYn — hE, fn = constant, 


In the explicit case, when Öz = 0, solving this equation means only a simple division 
of a constant expression with o. In the implicit case and a N-dimensional system 


19 


20 CHAPTER 2. LINEAR MULTISTEP FORMULAE 


function f we must, generally, solve a nonlinear system of N equations, so that it is 
necessary to use an iterative procedure for this special problem. 

The classical procedure for solving the implicit equation or the system of implicit 
equations is the predictor-corrector scheme. The predictor linear multistep formula, 
usually an explicit one, give a starting value for the iterative process for solving the 
implicit corrector formula. The approximation delivered by the predictor formula is 
improved applying the corrector formula. The standard iterative process for solving 
the implicit corrector equation is the method of simple iterations, in the nonstiff case, 
and the NEWTON’s method, in the stiff case. 

The formula (2.1) can be applied only if we know the values yo, ..., Yķ-1, referred to 
as the starting values. These are delivered by a procedure which is independent from 
the basic formula, a starting procedure. Naturally, the stepnumber of such a procedure 
is smallest that of the multistep formula. In almost all present numerical codes, the 
starting procedure includes some onestep methods. 

The starting procedure must satisfies two conditions: the boundless and the compa- 
tibility. The boundless refers to the existence of a supperior limit for the approximation 
values, the procedure outputs, an uniformely bound relative to a small variation inter- 
val for the stepsize. The starting procedure is compatible if the outputs values converge 
to a fixed value, when h — 0. These two condition ensure the existence and the unicity 
of the numerical solution. The linear multistep formula has an unique solution for a 
specific variation interval for the stepsize and for any bounded and compatible starting 
procedure. 


2.2 Consistency and order 


The consistency of a linear multistep formula makes sure that the numerical solution 
(the theoretical one) converges to the exact solution when h — 0. In this condition, 
we can find a stepsize h for which the approximation error is under some tolerance 
level. 


2.2.1 Definitions 


Let the linear operator associated with a linear multistep formula [86] 
k , d 
L:= lt “hö, (2.2) 


where H is the displacement operator defined as Hy(t) = y(t + h). 


Definition 2.1. L has the precision degree p, if L = 0 for any tf, j = 0,...,p, and p 
is the maximum integer that safisfies this property. 


The local error of a linear multistep formula is the difference 
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where y(t,) is the exact value of the solution at ta, and yy, is the approximate value 
of the solution at tn, produced applying the formula with x, = y(tn-i), i= 1,..., k. 
If the system function f is differentiable continous, then the local error can be appro- 
ximated by (MIRANKER[85]) 


yu), s (ou! - 2-0) — L(ylta — hk)), 


where p is a value between g(f,) and yn. Thus, we can justify the following definition. 


Definition 2.2. A given linear multistep formula has p order if the associated linear 
operator has a precision degree p. The linear multistep formula is consistent if its 
order satisfies p > 1. 


2.2.2 Algebraic conditions 


Let the expansion in TAYLOR’s series of the linear multistep operator 


L(y(t)) = do eh'y(t). 
2=0 
A precision degree p involves the relations c; = 0, 1 = 0,...,p. The values c; can be 


expressed in terms of a; and ĝi. 
The algebraic system, which is equivalent with the condition of p-order, is the 
following: 


Moreover, 


k 
Cprı = (p + (P+D! x Ak ett —(p + 1) ıı. 
0 


For a linear multistep with k steps, the inequation 


is referred to as the order barrier. Consequently, there is an only one implicit 
formula with k and of maximum order 2k and an only one explicit formula with k 
steps and maximum order 2k — 1. 


Remark. A p order linear multistep formula exactly solves an equation with a 
polynomial solution of lowest degree than p. 


The methods with k steps are generated by polynomials pairs (p, o) with real 


coefficients, where 
k k 
2435-0293 E. (2.4) 
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Suppose that (o, o) = 1. 
If E is the translation operator, specifical for the numerical process, Fy, = Yn41, 
then the linear formula with & steps can be rewritten in the form 


p(B) yn = ho(E) fx 
The consistency condition can be algebraic expressed in the relationships 
9(1) =0, 711) =o(1) 20. (2.5) 


The root € = 1 of the polynomial p is referred to as the principal root, and the others, 
as the external roots. 
The p order condition can be also set in the terms of the pair (p, o): 


ple”) — ha(e") = O(h?*), h—”0. 


Using the transformation e” +> x we get an equivalent equality 


Pe) (a) =Ol(e@-1)), L 


lng 
Note : 
z—l k z+1 
= (28) H) - Yoav. 
0 
z—l k z+1 k . 
s(z) =( 5 ) (5) de (2.6) 
j= 


The p order conditions can be expressed as the WIEDLUNG’s relationships: 


ak- ak- bk ak- bı 
1 
dı) 1 l —br, j odd . . 
5 7 bp ja + 30k-j-3 + xÖk-i-s free il „j € min(p, k) (2.7) 


- hue j even 


Suppose a; = b; = 0,1 < 0,2 > k. Therefore, 


The consistency condition can be rewritten: 
9(1) =a, =0,  pfl) = ol) = bi = ar /2 £0 (2.8) 


The difference between two methods with the same order can be overtake based 
on the error constant. The local error of a p order formula can be represented in the 


form 
ylin) — Yn = hP y OHD + O(t), (2.9) 
Definition 2.3. The value c* = —cp+1/0(1) is referred to as the error constant of the 


linear multistep formula (p, o) of p > 1 order. 
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The error constant can be easily computed (MIRANKER[86]): 


gə e = lool) 
rl (x — 1)Pr1 ` 


The p order accuracy supposes a p order of consistency and a small error constant 
(relative to the stepsize). 


2.3 Stability 


There are many examples of consistent numerical formulae which not converge to the 
exact solution when h — 0. The reason of this phenomenon is the accumulation of 
the calculus errors. The most refined is the integration division, the computational 
effort increases and it is possible that, at h — 0, the accumalated errors introduce 
some perturbations which are bigest than the exact solution values. 

Let consider an initial value problem and the system with some perturbated initial 
condition. A good numerical method for solving the initial value problem do not 
increase without limits the ratio between the actual difference of the two numerical 
solution and the initial difference (the perturbation of the initial condition), in any 
condition for the stepsize. This property is referred to as stability of the numerical 
method. 


2.3.1 Zero-stability 


The stability study of the linear multistep methods is based on the concept of root 
condition. 


Definition 2.4. A polynomial satisfies the root condition if all its roots are inside 
the complex closed unitary disk and those from the boundary are not roots of the 
derivative polynomial. 


Theorem 2.1. /85/ If a linear multistep method (p, o) is stable, then p(€) satisfies the 
root condition. 


The polynomial p is referred to as the stability polynomial. If p satisfies the root 
condition, we say that the formula is zero-stable. 
The reverse implication is true when the linear multistep formula is also consistent. 


Counterexample.[52] Let the consistent formula 
Yn + AYn—1 — DYn—2 — h(Af, + 2 fn—1)- 


The stability polynomial has the root —5. If we apply the formula to the differential 
equation y’ = —y, with y(0) = 1, we get the equation 


Yn + ACL + hyn — (5 — 2h) yn = 0. 
The solution of the this equation has the form 
Yn = axi (h) + bad (h), 
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where x;(h), i = 1,2 are the roots of the equations 


z? +4(1+h)x — (5 —2h) =0 


ie. zi(h) =1—h4+O(h?), a(k) = —5 — O(h). Note that |y,| — œ, n oo, 


The stability concept is extremely important for the convergence of the iterative 
process. The approximations produced by a consistent linear multistep formula con- 
verge to the exact solution when h — 0 iff the formula is zero-stable, i.e. 


Consistency \ 
Stability €— — > Convergence 


N, ( Root condition J) 7 


In the construction of efficient numerical codes, the DAHLQUIST 5 barriers are very 
important. These are the theorems which impose some restrictions on the stable 
formula class. The first DAHLQUIST’s barrier is the following theorem. 


Theorem 2.2. /53/ The order of a consistent and zero-stable linear multistep formula 
can not exceed the following values 


ET, ifk üs odd, 


+ 2, if k is even. 


Remarks: 


(1) The stability barrier is more restrictive than the order barrier. This fact is valid 
only for the linear multistep formula, not also for the onestep methods like the 
RUNGE-KUTTA’s process. 


(2) Moreover, HAIRER and WANNER [53] have proved that, p < k, if G,/a, < 0 
(particulary also for the explicit formulae). 


If all the secondary roots of p are inside of the complex unitary circle, the formula 
is strongly stable. A formula is strongly stable iff r(z) is an HURWITZ polynomial, i.e. 
it has all the roots inside the complex half-plane Rez > 0 (that implies a, > 0, j = 
0,...,k). 


2.3.2 Absolute stability. Stability domain. A-stability 


The absolute stable methods are particular cases of zero-stable methods. The idea is 
to request the decrease of the secondary numerical solution with the increase of the 
stepnumber. The study is made on the scalar test equation y’ = Ay, ReA < 0. The 
absolute stability is the minimum property that must be imposed on any numerical 
integration method. The behaviour of the method on the test equation is a prevision 
model of the behaviour in the case of some nonlinear systems. 
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Let the equation produced applying the method with constant stepsize h to the 


test equation: 
k 


S (a; — qb;)Yn-k+; =9, q Ah (2.10) 
0 
The characteristic equation can be get from this formula replacing the values yn-k+; 


with é: 
X aë -qY BiG =0, q=àh. (2.11) 
j=0 0 


The polynomial p — qo is referred to as the characteristic polynomial of the linear 
multistep formula. 


Definition 2.5. A method is absolute stable for a value q = hA if, for the value q, 
the characteristic polynomial satisfies the root condition. 


Theorem 2.3. //5/ A consistent linear multistep formula is absolute stable at q = 
hà € C iff all the solutions {yn} of the equation (2.11), produced applying the formula 
to the test equation yi = Ay with constant stepsize h, are bounded when n > ov. 


The proof of the direct implication is similar to the one of the implication stability 
— zero-stability. The reverse implication makes use of the consistency hypothesis. 


Counterexample. Let the same formula like in the last section and the equation 
yi = Ay, with y(0) = 1, ReA < 0. The solution of the difference equation has the form 
Yn = all — hà + O(h?d?)]" + B[-5 — O(h)]” where |hA| < 1. Thus |y,| — oo since the 
second values converge to infinity. On an other hand, also |1 —hA| > 1, Yk > 0. 

We can associate to any linear multistep formula a stability domain. 


Definition 2.6. The absolute stability domain of a linear multistep formula is the 
set 


A= {hA| applying the formula to the problem y’ = Ay, A € C, y(to) = Yo, 
with the constant stepsize h > 0, we get some approximate values 
Yn Which yn — 0, when n — co}. 


The stability domain is the set of the complex points for which the linear multistep 
formula is absolute stable. 
The boundary of the stability domain lies in some parts of the curve[53] 


— ple) 
q(0) 1 o(ef) 


The stability domains of the classical methods, like the ADAMS-BASHFORTH, 
ADAMS-MOULTON, NYSTRÖM or MILNE-SIMPSON’s formulae, are almost all 
bounded. The stability requirement leads to the condition that the values |hA,;|, where 
A; are the eigenvalues of the JACOBIAN’s matrix at the moment £f,, must be inside the 
stability domain. A large value |à;| implies a small value A. If the stability domain 
is unbounded, for example includes the complex left half-plane, then the condition 
hà; € A does not restricted the stepsize. 


0 € 10, 2). (2.12) 
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Definition 2.7. A formula is A-stable if the absolute stability domain includes the 
complex left half-plane, i.e. those values z for which Re(z) < 0. A method is A-stable 
if it is based on an A-stable formula and, moreover, the implicit equation produced 
applying the formula at the test equation is exactly solved. 


Remark. The method notion is more general than that of a formula, scheme or 
process. Especially in the stiff case, a method must includes also a procedure for 
solving the implicit equations. We can also include in a method the step variation 
techniques and the order variation tehniques, the estimating error algorithms or the 
algorithm for stiff detection. In the following sections, we refer to a method only as a 
formula, scheme or process in an implementation with a constant stepsize. 


Let 
C_ = {z € Cİ Refz) < 0}. (2.13) 


Thus, the condition of A-stability can be expressed by the inclusion 


and that of zero-stability by 
Ü € A. 


For an onestep method, the A-stability condition is derived from a property of 
the stability function R(q) = Yn+1/Yn. An onestep method is A-stable if the stability 
function satisfies |R(z)| <1, VRe(z) < 0.. 

Also, if the linear multistep formula has an unbounded stability domain, it is 
possible to have an undesirable property: 


lYnti/Yn| > 1, hà => —oo. 


For example, the trapezoidal rule, which is A-stable, has this property. Thus, the fastly 
decreasing components of the exact solution are represented in the numerical solution 
by some slowly decreasing components, even by some oscilatory ones, if Yn41/Yn < 


0, Vn > 0. 
Definition 2.8. 


(i) A linear multistep formula is stable at infinity if limpy+—co |Yn41/Yn| = 0, where 
Yn is the numerical solution of the scalar test equation. 


(ii) The formula is L-stable or stiffly A-stable if it is A-stable and it is stable at infinity. 


(iii) A linear multistep formula is strongly A-stable if it is A-stable and there is a 
value w > 0 such that 
sup |R(hA)| <1 


h\<—w 
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2.3.3 Sufficient conditions for A-stability 


A practical method to decide that a formula is A-stable consists in applying the fol- 
lowing lemma proposed by MIRANKER [85]. 


Lemma 2.1. If £,? = 1,...,k are the roots of the characteristic equation, then the 
following statements are equivalent: 


(i) the method is A-stable; 
(ii) if Req <0, then |€&| <1, i=1,...,k; 
(it) YE: |El > 1, Req(£) = Re[p(€)/o(€)] > 0 holds for |é| > 1. 


Note W the exterior of the closed unitary complex disk. Using the lemma 2.1 we 
can prove the following propositionİS61. 


Proposition 2.1. 


(i) The linear multistep method is A-stable only if the roots o, of o(6) safisfy lcil < 
l, i=l... k; 


(ii) The linear multistep formula is A-stable iff Re[p(€)/o(€)] > 0, v€ € W. 


(iii) If a root of o(6) lies on the unitary circle and is not a simple root, then the linear 
multistep formula is not A-stable. 


Applications: 


(1) The trapezoidal rule is A-stable since p(€) = € — 1, o(€) = 1(€ + 1), Req(£) = 
(161? — 1)/1£ + 1İ”and thus Re q($) > 0, for all |€| > 1, and the root of o from the 


unitary circle is a singular one. 


(2) The EULER’s implicit rule is A-stable since p(€é) = é — 1, a(€) = £, Re q($) = 
(Aİ — Ref)/İ£ + 1]? > 0, Vig] > 1; 


(3) An example of A-stable and consistent method for arbitrary stepnumber k > 1 is 
the following: 


kh 


Since p(€) = €” — 1, o(€) = Eek +1), the roots of o(€) are separated ones and 
lie on the unitary circle, and Regq(€) = k(qP" — 1/16 --117 5 0, Vİ/İ > 1. 


Proposition 2.2. /86/ If 
(i) the roots o, of the o polynomial satisfy |o;| < 1, i= 1,...,k, and 
(ii) u(€) = Reg(€) > 0 for € on the unitary circle, 


then the linear multistep formula is A-stable. 
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Applications: 


(1) Let the class of onestep methods which depend on the parameter a (referred to 

as a-rules) 
Ynt1 — Yn = hid m a) fra + a fal. 

We apply the last proposition. The unique root of the o polynomial is o) = 
—a /(1 — a) with |o,| € 1 iff a < 1/2. An other hand, u(e) = |o(e’’)|-? P(e’’), 
where P(e”) = (1 — 2a)(1 — cos 8). Then P(e") > 0 iff a x 1/2. Consequently, 
the method is A-stable if a < 1/2. Note that we have only an implication. A 
counterexample is that of the trapezoidal rule which is an A-stable formula with 


a= 1/2. 
(2) The two step formulae of order two can be represented in the form 
(1 — c)ya + 2cyn-1 + (1 — c)yn-2 = hibf, + (2 — a — b) fn-1 + a fn-2], 
where c= a — b. From Proposition 2.2, the inequalities 
b—asÜ, -l+a+6>0, 
are sufficient conditions for A-stability. 
Corollary 1. /f 
(i) Re(s;) < 0,2 =1,...,k, where s; are the roots of the polynomial s(z), 
(ti) Re(r(z)/s(z)) 2 0, Vz: Re(z) =0, 


then the linear multistep formula is A-stable. 


The request condition (7) can be treated with the ROUTH-HURVVTTZ”s theory. More 
details can be found in [52], [45] and [96]. 


Note 
k 


X ey” = Re [r(iy)s(—iy)]- 


0 
A sufficient condition for (ii) is e; > 0. The coefficients e; can be computed from the 
following relationships built up in the conditions a, = b; = 0, ¿ < 0, ¿i >k: 


Ee = NO (=1) H ajbai-j, ? =0,...,k 
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2.3.4 Necessary conditions for A-stability 


The class of A-stable linear multistep formulae is not very large, since such a formula 
must be implicit and with an order of accuracy lowest than three. 


Theorem 2.4. (ENRIGHT/38/) A linear multistep formula, with a stability domain 
which includes the negative real half-axis, must be an implicit one. 
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Corollary 1. An ezplicit lincar multistep formula cannot be A-stable. 


This result shows that it is necessary to solve at each step an implicit equation or 
a system of implicit equations (a great computational overhead), which can be also a 
source for suplimentary errors. 

The following theorem, due to DAHLQUIST, gives a full characterization of the 
restrictions impose by the A-stability condition. It is referred to as the second 
DAHLQUIST’s barrier. 


Theorem 2.5. (DAHLQUIST/34]) The order of an A-stable linear multistep formula 
can not exceed two. The trapezoidal rule is the two order A-stable formula with the 
lowest error constant ce” = 1/12. 


Corollary 1. For a linear multistep formula 


The are many strategies to overtake the second DAHLQUIST’s barrier. First of all, 
we can apply some RUNGE-KUTTA’s process. These must also be implicit when we 
impose the condition of A-stability, and at each step we must solve a system of q x N 
implicit equations, where q is the number of intermediary stages (compared with only 
N equations in the case of a linear multistep formula). The advantage of using such 
process is that of the identity between the barrier order and the barrier imposed by the 
A-stability condition on the accuracy order. The onestep method are not the subject 
of this study. More details can be found in [53], [16] or [96]. 

The principal ways to break the DAHLQUIST’s second barrier are: 


. weaking the stability condition: stiff stability, A(a) stability and so on; 
2. using the implicit RUNGE-KUTTA’s process; 
. using the high derivatives of the solution: 
a) multiderivative multistep formulae; 
b) multiderivative RUNGE-KUTTA’s process; 
. using some nonlinear interpolators: 


a) nonlinear multistep methods; 
b) rational RUNGE-KUTTA’s process; 


. adding some new stages, extradivision points or future points: 
a) cyclic linear multistep methods and composed methods; 
b) hybrid methods (extradivision points) and multistep 
RUNGE-KUTTA methods; 
c) extended methods (future points); 
. reducing the approximation error: 
a) exponential fitting; 
b) extrapolation methods; 
. exploiting the special structure of the system: stiff separability. 


(2.14) 


These ways for breaking the barrier are presented in the following sections. 
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2.3.5 A(a)-stability 


Small variations of the system JACOBIAN’s matrix produce small variations of its 
eigenvalues A;. By such a variation, the values hA;(t) cover only a small part of the 
complex left half-plane. The covered region can be delimited by two lines with the 
slopes less than one. For the problems for which the JACOBIAN’s matrix eigenvalues 
are not included in the imaginary axis, it is justified the request of absolute stability 
on a smallest region than the entire left complex half-plane. 


Definition 2.9. A linear multistep formula is A(a)-stable, with 0 < a < 7/2, if all 
the solutions of the difference equation, produced applying the formula to the scalar 
test equation y’ = Ay, converge to zero when n — oo, for any fixed stepsize h > 0, 
and for any A Æ 0 so that 


hà € Cy = (ql larg(—q)) € o, q £ 0} 


holds. 


The angle a is referred to as the stability angle of the method. 
Formally, A(a)-stability means 


C. CA. 


The algebraic conditions can be deduced similarly to the case of the A-stability 
property: 


e a linear multistep formula is A(a)-stable if, for any q € Ca, all the roots of the 
characteristic equation satisfy |&| < 1, 7=1,...,k. 


ə a linear multistep formula is A(a)-stable if all the roots s, of the polynomial s(z) 
satisfy Res; < 0,2 =1,...,k; 


e a linear multistep formula is A(a)-stable iff r(z)/s(z) lies in the exterior of the 


set Ca, Vz: Rez > 0; 
The most important properties of the above notion are the followings: 
1. A(a)-stability = A(@)-stability for 0 < 8 < a; 
2. A(z /2)-stability = A-stability; 


3. an explicit formula can not be A(a)-stable, since C, includes the negative real 
half-axis; 


4, it is possible to overtake the second DAHLQUIST’s barrier, since are some 
A(a)-stable formula for 0 < a € 7/2 and k = p = 3 or k = p = 4. 


Theorem 2.6. (GRIGORIEFF, SCHROLL/ÖöT/) Leta < m/2. Then for any k a positive 
integer number there is a A(a)-stable linear multistep formula with k steps and of 
order p = k. 
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The limit case is a = 0. 
Remarks: 


(a) An explicit formula can not be A(0)-stable, since the stability domain of a such 
formula must includes the negative real half-axis; 


(b) The trapezoidal rule is the only one A(0)-stable formula with p > k + 1. 
For a linear multistep A(0)-stable formula a; > 0 (or a; < 0) and 9, > 0 (or 
B; < 0) hold for all 7 = 0,..., k (LINIGER[81]). 
2.3.6 A%—, 40—, A,.—stability 
Building up these notions, the idea is to weak the condition of A-stability. 
Definition 2.10. Let a linear formula with k steps. 


(i) The formula is A®°-stable if, applying it to the test equation with A a negative real 
number, we get some approximations yn with the property yn > 0 when n > oo, 
for any constant stepsize h (NEVANLINNA İ90)). 


(ii) The formula is Ao-sfağle if the characteristic polynomial satisfies the root condi- 
tion for any negative real value q (CRYERİ32İ). 


(iii) The formula is A..-stable if s(z) is an HURWITZ polynomial (LINIGER[81]). 
The most important properties are the followings: 

(1) A®-stability = (—œ,0] c A. 

(2) A(0)-stability => Ao-stability => A®-stability. 


(3) The formula is Ap-stable iff all the roots of the polynomial r(z) — qs(z) are inside 
the set C_ for all q € (—oo, 0) (CRYERİ33İ). 


(4) If the formula is Ap-stable, then a;b; > 0, j =0,...,k. 


(5) The formula is Ao-stable iff Vz € C4, r(z)/s(z) is regular and positive, and for a 
z value with z = Imz, r(z)s(z) > 0. 


(6) If the formula is A..-stable, then b; > 0, 7 = 1,...,k and the order p < k 
(LINIGER[81]). 


(7) The formula is Ap-stable if it is Az-stable, ao > 0, and the polynomial 


k-1 
1 , 
uni (iy) (-üylz So fi, Ve R, ny” xü, 
0 


has only real positive roots. The last condition holds if f; S 0, 7 = 0,...,k— 1. 
(8) The formula is A(0)-stable iff (JELTSCH[66]): 
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(i) it is Ao-stable: 
(ii) all the roots €; of the polynomial o(€) with |&| = 1 satisfy o’(€;) 4 0 and 
Re (p(&)/(&:o"(&i)) > 0; 
(Hi) all the roots £ of the polynomial p(€) with |€;| = 1 satisfy 
Re (o(&;)/(&p'(&;)) > 0. 
Application. The methods with k steps, the order k, and the polynomial s(x) = 
(a + diy (the CRYER’s formulae), are Ao-stable iff d > Qk+1 | 


Countererample. The formula 


h 
Yn — Yn-1 = ri: + 2fr-1 + fn-2) 


is Ao-stable, but not A(0)-stable. 


2.3.7 Stiff stability 


Let a linear multistep formula and €;, the roots of the characteristic polynomial. Then 
S = {à € CJE) a. (2.15) 


is referred to as the relative stability domain, and also as the star region. 


Definition 2.11. (GEAR[45]) A linear multistep formula is stiffly-stable with respect 
to the D, a, 0 parameters, if 


Hu U hə CA, R3 CS (2.16) 
where 
Ry = (A € Cİ ReA < —D}, Rp = {AE C| — D < Reà < —a, |ImA| < 9}, 


R3 = {A € Cİ |ReA| < a, |ImA| < 9). 


The reason of this notion is the following. For the test equation the exact solution 
increases at one step with €”. If BA = u+ iv, the amplification term is e”. If 
u € —D < 0, the solution component decreases very fast. Since the solution is 
very small, we are not interesed in accuracy of that, only in its stability. In the 
neighbourhood of the origin, the accuracy becames important. The relative stability 
is a necessary condition for accuracy. If |v| > 6, at least 9/(27) complete cycles are 
performed at each step. These variation is very hard to be followed if @ is large, and 
we can use a reasonable stepsize (relative to the integration interval length). 

The classical examples of stiffly stable formulae are the backward difference formu- 
lae of maximum order (GEAR’s formulae). 

The stiff-stability notion is very important in breaking the second DAHLQUIST’s 
barrier. For any arbitrary stepnumber & we can build a stiffly stable linear formula 
with k steps. 
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Theorem 2.7. /CRYER/32/) For any value D > 0 and any positive integer k there is 
a linear formula with k steps, stiffly stable, and of k order, with 


Rı CA, 


where A is the stability domain of the formula. 


The proof is based on the fact that for any D we can get a value d such that the 
relationship from the theorem holds for the CRYER’s method generated by o(€) = 
(E+ diy and the condition of maximum order. 

The relationships between the stiff-stability notion and the above mentionated 
notion can be expressed in the following diagram: 

A-stability => stiff stability = A(arctg(@/D))-stability > 
= A(0)-stability => Ao-stability 

Coutererample. The method generated by the pair of polynomials 


pl) = — €?/2 — 1/2, ($) = 36/2 — © + 3€/2 
is A(0)-stable but not stiffly stable. 


A reverse implication supposes some supplimentary conditions: 


Proposition 2.3. (JELTSCH/66/) A linear multistep formula is stiffly stable iff 
(i) it is Ag-stable; 
(ii) the roots of the polynomial p(€)/(€ — 1) are inside the unitary circle; 


(iti) all roots &; of the polynomial o(€) for which |€;| = 1 satisfy o'(€;) 4 0 and 
P(E) /(io'(&i)) € Ry. 


Note that COOKE[30] has established some sufficient conditions for which a con- 
sistent linear multistep methods is stiffly stable. 


2.4 Convergence 


The consistency condition supposes that the numerical process exactly models the so- 
lution of the initial value problem when the stepsize h — 0. One principal supposition 
consists in the fact that the calculus is carried without any perturbations. The more 
general notion of convergence makes reference to the same phenomenon, but allows 
the intervention of some small perturbations. 

Any numerical method must satisfy the condition of convergence to demonstrate its 
efficiency in the integration of an initial value problem, since the convergence supposes 
a control of the error produced by the numerical process. 

The global error (accumulated error) of a numerical method is 


En = (tn) — Un, (2.17) 


where y(t,,) is the exact solution, and oy is the approximate solution at tn, produced by 
the numerical code. The boundary of this error leads us to the notion of convergence. 
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Definition 2.12. A method is convergent if 
lim max İle, İİ = 0 (2.18) 
h>0 n 


holds for all the initial value problems with the standard conditions for the solution 
unicity. 


2.4.1 Global error of an onestep method 


Let the onestep method 
Yntl = Yn + hö(h, Yn, Yn+1)- 


Note the global discretization error 


Enti = Y(tnti) — Yni, (2.19) 


produced by the accumulation of the local discretization errors, 


Tn = (tası) — Una” 


where 
Yaa = Y(tn) + hO(h, yahu). 


The onestep method is of p order of consistency, if p is the greatest integer which 
satisfies 


İsa = O(P, h—0. 


If the method is stable, then there is a constant value K for which 


naa Z Yill < AUG, — Yall, Yn = y(n) 


holds, and R 
lenll < Kfe + nti || < Klel +2, 


where p is the superior bound of the local discretization error. Thus, the global 
discretization error €,4, satisfies the inequality 


[Enli KH ol) +S Ki. (2.20) 


i=0 


If K < 1, he. the method is contractive, there is not a significant propagation 


of the global error, only an accumulation of the local error, and, if €, = 0, then 
enll S 7/01 — K). 

If K = 14+ Coh, h € (0, hol, where ho, Co are independent of the initial value 
problem, then 


_ qz , (+ Coh)”ı—1 
Eel SPY OU + cay - pA < 


i=0 
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pho tnat, Co <0 
— m. — 1), Co >0 ? h € (0, ho] 


when İleoli = 0. 

If there is not a constant value Co so that K = 1 + Coh, then the method is not 
stable. 

If we approximate the exact solution with the perturbated numerical values v,,, the 
propagation of these perturbations can be bounded in the same manner like in the 
propagation of the global discretization error. 

Note the global numerical error 


Čn+1 = azi — Unt; (2.21) 


produced by the accumulation of the local numerical errors 


la41 = Uni — n41, 


where 
Unt = Un + hö(h, Ün, Yn41)- 
We get 


lëna < Klè ll 2> lhal leli 0223. 


i=0 
Hence, the global error 
Enti = Yltn41) — Undi = Engi + Enti, (2.23) 
satisfies a 
lensill < KH AlE + lleoll) + (B+ A) 2 A’. (2.24) 


2-0 


In practice, we suppose that p < P. 


2.4.2 Convergence of the linear multistep formula 


Note y,(t) the interpolation function of the approximate values, produced applying a 
linear multistep formula at an integration interval division with constant stepsize h. 


Definition 2.13. A linear multistep formula is convergent if 


lut) =H +0, h= 0, £ € [to to + Tİ, 


holds for any initial value problem with the standard conditions for the unicity of the 
exact solution, independent from the starting values, which only satisfy 


lylto + ih) — y(to +ih)| 20, bh 0, i=0,...,k—1. 


The algebraic conditions for convergence are the zero-stability and the consistency. 
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Theorem 2.8. /52/ If a linear multistep formula is convergent, then 
(i) is zero-stable; 
(ii) is consistent. 


The reverse implication is based on the transformation of a linear multistep formula 
into an onestep method. Let Y; = (yj-1, yi-2,---, Vin), 1 > 0, 


—ap-1/A, —Ap-2/Op ... —Oo/dik 
1 0 2. 0 
A= .. .. : 
0 1 0 


and the function Y, uniquely defined by the relationship, 


Br va arj 
ti— s 00— ti, Aw — 1— R 
ə Bi F( q? Mi j) Tü ot ə On Yi-j 


j=0 


Then the linear multistep formula can be rewritten in the form of an onestep formula 
[53]: 


where A © / is the KRONECKER’s tensorial product, and 
O(t, Yi, h) = (ei © DD (t, Yh), ci = (1,0,...,0)". 
Theorem 2.9. /53/ The zero-stability and the consistency imply the convergence of 


the linear multistep formula. 


Moreover, the global error can be estimated when the formula is A-stable or only 
A(a)-stable. For a linear system 


yi = Ay + g(t). 


the following result holds. 


Proposition 2.4. /53] If the linear multistep formula has p order of accuracy, is also 
A(a)-stable, stable at infinity, and the system matrix A has the property that there is 
a matrix T so that TOAT = diag (à1,..., Ax), where |arg(—A;)| <a, i= 1,..., N, 
then there is a constant value M, which depends only on the formula, so that for any 
stepsize h > Ü, the global error satisfies 


lya) — g))ll s MIPIT (mə Gü) “gil ef m5522 (2.25) 


More generally is the following result concerning the nonlinear systems. 
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Proposition 2.5. /53/ If the method is strongly A-stable, of p order of accuracy and 
the initial value problem satisfies the disipativity condition 


Re € f(t,u) — f(t,v),u—v >< ullu — olf, 
then there is a constant value Co > 0 so that 


lyn = yita)ll s C(max llus — yt) + A max [Fn yi) = y's) + Mh”, (2.26) 


O<j<k 


holds if hu < Co, where C and M depend on the method, and, for > 0, also on 
t, —to. Moreover, M depends on the exact solution derivatives of a given order. 


2.4.3 Convergence of the discretization methods 


CHARTRES and STEPLEMAN[27] present a theory of convergence, consistency and 
stability in which the definitions are more general than the above mentionated ones. 


Definition 2.14. A discretization method is a set (F,{Fi,},{E,}, H, X,Y) where 
Fix >Y, Fa: X x FE, >Y for each h € H, 0 € Ek, in H we have defined the 


convergence h — 0, but 0 ¢ H, and Y is a linear normed space. 


In this definition we have note: 
ə X: the date space, the set of pairs (yo, f); 


ə Y: the function space which includes the exact solution and the interpolation 
function of the dicrete values; 


ə F: X — Y maps each pairs (yo, f) € X into the exact solution y € Y; 


ə Fp: the set of the N, components of the perturbations produced at each step 
i = 0,..., Np of stepsize h (practical, is a finite dimensional space); 


ə Fi: X x E, o Y maps each pairs of an element from X and a specific pertur- 
bation into the interpolation function of the numerical discrete values; 


e H: the set of the admissible stepsizes. 


Suppose, for each method and each h € H, the existence of a function öy, : Ep > R, 
with ¢,(0) = 0 for any h € H. Then by e, — 0, we mean (en) > 0, and by 
(h,e,) 20, h +0 and dalen) — 0. 


Definition 2.15. 


(i) A discretization method is convergent at x if liM zo Fh(£, en) = F(x). The 
order of convergence is p if İF) (a, ea) — F(a)ll = O(h’), for da (ea) = O(h’). 


(ii) A discretization method is consistent at x if there is a specifical perturbation ef 
for which h € H such that limao pa (e$) = 0, and limp4o Fa (ə, ek) = F(x). The 
order of consistency is p if İF (a, e$) — F(x)|| = O(h?) for o, (ef) = Ohe). 
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(iii) A discretization method is stable at x if lim(,..,)+0 İF (ei en) — F(x, 0)|| = 0. 


Theorem 2.10. /27/ A discretization method is convergent iff it is consistent and 
stable. 


Application. The EULER’s explicit rule applied to an initial value problem on the 
integration interval [0,1] is a discretization method (F,{F,},{ En}, H, X,Y) with 


X= {luo AF : 10, 1) x R => R, for which y — F(t,y), y(0) = Yo; 


has an unique solution y}, 


Y= Cl, 1], F(f, Yo) = y, En = RYH, 


N, = L, H = {1,1/2,1/3,...},  Fil(f,yo),e) =u 
where 
u = İy(uo, u1, ..., UN, ) 
uo = Yo + €o, Un = Uni + hflta-1,Un-1)+ en, 2 =1,2,...,Np (2.27) 
e = (eo, €r,.... ex, ) € En 
If u* = I,(ug,---,uy,) is the numerical solution associated to en = 0,n = 0,..., Nn 


and f satisfies the LIPSCHITZ’s condition in the dependent variable, then 
[uo — up|] < eo [um = wp] < Q + AL) en = uzali + İle: İİ, 
where L is the LIPSCHITZ’s constant, and, therefore, 
lun ue See lel, n= 0)... Nİ. 
r=0 


Hence, the most weakest topology in Ep, for which the method is stable, is based on 
Pale) = SA le,|. Note TE, the local discretization error in £,. Then 


y(tn) = y(tazı) + hf(azı, y(tn-1)) + TEn 


holds and the characteristic perturbation ef is the vector (0, TE1,...,TEn,). In the 
proposed topology (e$) = Yu T E,. The method is consistent if ġa (ek) — 0, when 
h — 0. Since TE, = OÖ(h?), if f is differential continous and has a bounded derivative, 
the method is convergent. 


2.5 Linear multistep formulae for stiff problems 


There are three ways to build-up a linear multistep formula: 


e by numerical integration; for example, we can mention the ADAMS-BASHFORTH, 
ADAMS-MOULTON, NYSTRÖM, MILNE-SIMPSON or NEWTON-COTES formulae; 
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e by numerical diferentiation; for example, we can mention the backward differen- 
tiation formulae; 


e by combination of the coefficients of some known formulae. 


A bounded stability domain is a characteristical property for the methods of the 
above mentionated first way. Consequently, these formulae are not suitable for the 
integration of stiff problems. Most adequate are the following methods: 


e the backward differentiation formulae (short noted BDF). For k < 6, the methods 
with & steps and maximum order k, known as GEAR’s formulae, are stiffly stable; 


e the methods of ADAMS’s type. Note C} the ADAMS-MOULTON”s method of q+ 1 
order with q steps, and the generating polynomials 


q 1 
pale) =E, ale) =X yle), = o f Cidu, 
j=0 0 


Let 6,,...,5, some real values. We can associated to those values a, method of 
q+ 1 order with q + k steps (RODABAUGH, THOMPSON{118]) 


k k 
C (bi, pb) = So Com + (: -> 2 Cy; (2.28) 
m=l1 m=1 


e the a-methods with the general formula 


k 
Yn — WYn—-1 + (a — 1)Yn-2 = 33577 Ü co o 2. (2.29) 


0 


For a = 1 we get the ADAMS’s methods. ROCKSVVOLDİ117İ proves that there 
are such methods which are Ao-stable, for the orders p = 1,...,4; 


e the CRYER S methods are defined by the polynomial o(e) = (¢+d)*, d € (-L1) 
and the condition of maximum order of accuracy (CRYERİ321). JELTSCH[67| 
proves that the CRYER”s method with k steps is A(0)-stable and stiffly stable if 
—1 < d x —1 +2/(1 +2**'), where k = 1,...,7; 


e the GRIGORIEFF-SCHROLL’s methods [51] are computed using the polynomial 
o(e) = (e+ de +c) and the condition of maximum order. The numerical 
tests have proved that, for some values d and c, the methods are stiffly stable 
for k < 12 steps. 


2.6 Backward differentiation formulae (BDF) 


The BDF are the first formulae used in the integration of the stiff problems. 


40 CHAPTER 2. LINEAR MULTISTEP FORMULAE 


2.6.1 Formula 


A backward differentiation formula with k steps, BDF;,, can be represented by the 
iterative process 


k 
3377 


j=0 (2.30) 
where (az + 8?)ao Æ 0. 


The formula of maximum order is the result of a numerical differentiation proce- 
dure. The basic idea is to determine an interpolation polynomial for which 


P(trti-k) = Ynti-k; i= 0, ... ,k, P' (ta) = F (tas Yn) 


holds. The polynomial is built-up using the backward differences. Note Į the back- 
ward differentation operator. Thus, the GEAR’s formula with & steps is the following 
process 


Roa 
j=l 7 


2.6.2 Particular cases 


e The GEAR’s formulae [45] or BDF with maximum order (equal to the stepnum- 
ber), are stiffly stable for k < 6 steps. 


e The VARAH’s formulae [120] reduce with an unity the method order, in the 
favour of a greatest stability angle than that of GEAR’s formulae. 


ə The CASH’s extended BDF [22] suppose an implementation in a predictor- 
corrector scheme: the predictor formula is one of the GEAR’s methods, and 
the corrector formula is an extended one: 


k k 
(P) 2./ üilsi-k — hBy fas (C) ə OGUn+tj—k = hök. fn + hGrarfn4i- (2.32) 


j=0 


The scheme has k + 1 order of accuracy. The purpose of using this scheme is to 
improve the maximum order of the stiffly stable formulae, and to maximize the 
stability angle. One step consists in the following stages: 


1. compute y, by (P); 
2. compute 9,4, by (P); 


3. evaluate fi.443 


4. compute yn by (C). 
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2.6.3 Similar methods 


e The split linear multistep methods have been proposed by CASH [24]. A split 
formula with k steps has the form 


k k-1 
2. ings — h 2. BiFnt ik = ROL (tT) + (Be — A) F(tns Yn) (2.33) 


1-0 1-0 


where 7, is a value generated by a predictor formula of the form 


k k 
(P) X angie mh) B; Farin 
j=0 


j=0 
CASH has studied the effect on the stability of splitting the BDFs. 


e The DILL’s formulae [36] lie in the iterative process 


k 
N ajYnyj=k = h( Br fn + Bn-ufn-u), Pk = l, u€ {1,..., k} (2.34) 
j=0 
The purpose of building these formula is to increase the order of stiffly stable 
formulae. 


2.6.4 Stability 


A suggestive proof for the stifEstability of GEAR’s formulae is the study of the stability 
domain plots. The boundary of the stability domain for the k-step formula lies on the 
curve 


Thus, for k = 1, 
Az (sll:-1121) 
and for k = 2, 
OA = {z|Rez 3/2 — 2 cos 6 + 1/2 cos 26, 6 € [0,27)}. 


hold. The stability domains for k < 7 are plotted in Figure 2.1. 

The coefficients and the parameters of stiff-stability are given in Table 2.1. 

Note that the first GEAR’s formula is the implicit EULER’s rule. 

HAIRER and WANNER [55] have proved that the GEAR’s formulae are unstable for 
k > 7. 

The establishment of the stepsize variation interval is a very important problem. 
The stiff stability property can be lose applying an unadequate stepsize variation 
technique (GEAR et al.[46] [47]). 

The VARAH’s formulae [120] increase the order of the stiffly stable formulae with 
two units. The formulae are A(a)-stable for k < 12 steps and p < 8 order. The 
stability angles are mentionated in Table 2.2. 

The CASH’s extended schemes have two qualities: 
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Figure 2.1: Stability domains of the GEAR’s formulae 


Table 2.1: The GEAR’s formula coefficients,the stability angle and the stiff parameter 


1. the maximum order of a stiffly stable formula is 9, greatest than 6, the maximum 
order for GEAR’s formulae (see Table 2.3); 


2. it proves the possibility to break the DAHLQUIST’s second barrier, since the 
scheme of p = 3 and p = 4 order are A-stable. 


An other example of breaking the second barrier is the following CASH”s split 
backward differentiation formula: 


_ 18 5 9 2 3 6 — 
(P) Jn- (= + 20) Yn-1 + (= + 2 Yn-2 — (= + z0) Yn-3 = (= — 0) hf,. 


18 9 2 6 — 


The scheme has p = 3 order and it is L-stable for 9 € (—1.1, —0.29). In the same 
manner, we can built-up some split backward differentiation formulae with k > 3 
steps, in the idea to improve the stability angle of the GEAR’s formulae. CASH proves 
(by plotting the stability domains boudaries) that: 
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Table 2.2: The stability angles for the VARAH’s formulae 


T 8 9 10 11 12 
5 6 7 7 8 8 
a 17255 29 39 2 30 


Table 2.3: The stability angle for the CASH’s extended BDF schemes 


1 2 3 4 5 6 7 8 
BE 3 4 5 6 7 8 9 
o 190 90 90 87.61 80.21 67.73 48.82 19.98 


- the 4-step formula of 4/h order is L(a)-stable with a > 89.999° for a fixed 6 value, 
- the 5-step formula of 5 order is L(a)-stable with a > 85° for a fixed 9 value, 
- the 6-step formula of 6/h order is L(a)-stable with a > 71° for a fixed 6 value. 


An other example is the split scheme based on the trapezoidal rule: 


(P) Yn — Yn-1 = Ohf (ta, Un) + (1 ~~ Ah f (ti, Yn-1) 


(C) Yn — Yn-1 = Ohf (tn, Yn) + Sh F(t Un-ı) +h G — 0) Flin Yn) 


The scheme is stable at infinite if 9 = 1 + 1/V2, A-stable if 0 > 1/4, and we get the 
minimum error when 6 = 1 — 1/V2. 

The stiff stability of the DILL’s formulae have been also proved in [36] plotting the 
stability domains. 

The stability results concerning the above mentionated methods are concentrated 


in the following lemma. 


Lemma 2.2. Let k the stepnumber of a backward differentiation formula. 


(i) The GEAR 5 formulae are A-stable for k < 2 and stiffly stable for k < 6. The 
mazimum order for a stiffly stable formula is 6 (GEAR[{5]). 


(ii) The VARAH’s formulae are stiffly stable for k < 12. The maximum order for a 
stiffly stable formula is 8 (VARAH/120)). 


(iti) The CASH 5 extended BDF schemes are A-stable for k = 1,2 and stiffly stable 
fork <8. The mazimum order for a stiffly stable formula is 9 (CASH/22)). 


(iv) The CASH s split BDF schemes are A-stable for k = 3 and some 0 values, and 
stiffly stable fork <6 and some 6 values. The maximum order for a stiffly stable 
formula is 6 (CASH[24]). 
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(v) The DILL’s class of formulae of p = k — 1 order includes some methods of 7 and 
8 order which are stiffly stable. The maximum order for a stiffly stable formula 


is 8 (DILL/36/). 


Chapter 3 


GENERALIZED LINEAR 
MULTISTEP FORMULAE 


Building up some linear multistep formula of high order for the stiff problem integration 
we can follow two directions: 


1. using the high derivatives of the exact solution; 


2. adding some supplimentary stages, extradivision points or future points, steping 
in a large field of new methods like: 


(a) block methods; 
(b) hybrid methods; 
(c) extended methods. 
In the last section we have mentioned an example of an extended method. 


The low order linear multistep formula proposed in the last chapter can be im- 
proved, from the point of view of the error, applying one of the following methods: 


e exponential fitting; 
e extrapolation. 


In the following sections we discuss these ways. 


3.1 Second derivative linear multistep formulae 


The linear multistep formula uses only the derivative of the exact solution y’(t) = 
f(t, y(t)). A second derivative multistep formula uses also the second derivative of the 
exact solution y”(t) = d( f(t, y(t))/dt. 

The motivation for using such a formula is due to the fact that a linear multistep 
method for integrating the stiff problems must include a scheme for solving the implicit 
equation systems. The usual predictor-corector scheme can not be applied since the 
algebraic condition for convergence (||h(Of/Oy)|| to be a small value) is contradictory 
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with the stiff character of the problem. The iteration schemes for solving the implicit 
equations are based on some modifications of the standard NEWTON’s iteratives, which 
are exactly for linear systems. These methods avoid the difficulty of convergence, but, 
unfortunately, they ask the evaluation of the JACOBIAN’s matrix. Therefore, we can 
use this matrix explicitly in the linear formula. 


3.1.1 Formula 
Note g(t) := y” (t) = KƏ f föt) + (Of /Oy) f\(t, y@)). The iterative process, which repre- 


sents the formula, is 


k k k 
ə Ok-iYn-i — h ə Bri fai +h? ə Vk-i9n-i- 
0 0 0 


3.1.2 Particular cases 


e The ENRIGHT’s formulae are selected based on three requests: 


1. stability at infinity: a single nonzero coefficient y;, respectivelly yx; 


2. some reasonable stability properties in a neighbourhood of the origin of 
the complex plane: the polynomial p is the same like that of the ADAMS’s 
formulae. The iterative process is 


k 


i=0 


3. maximum order: the coefficients are unique defined (ENRIGHT [38]), 


j=l J i=j i=0 
where 


Ci = 


il (s = D’s(s + 1)s:2)(610-2) | 


il 
e The second derivative backward differentiation formulae (SDBDF) [53] have the 


coefficients derived from the condition of maximum order: 


: (> z) y’ Yn = (> z) hfa — ə. (3.3) 
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3.1.3 Similar methods 


e The CASH s second derivative extended schemes [21] consists in ENRIGHT’s me- 
thod as the predictor formula and an extended method as the corrector formula: 


k 
) ə Aj Yn+j— (ax, fanti- bth? So aida ky 


0 


k+1 


| Lom k = m... ath its ke (3.4) 


e A modification of the above mentionated extended formulae is also presented by 
CASH [18] as a composed second derivative multistep scheme, for autonomous 
initial value problems with f = Jf, where J = (Of/Oy) is the JACOBIAN’s 
matrix: 


k k k 
(P) Ym — Yn =h STi frari-e +hICD | Biynti-k +h >) Tifosi): 


j=0 0 0 
k+1 
(C) Ym — Yni = hd ay tng je + hd 7 755 b). (3:5) 
0 0 0 


The idea of this modification consists in the fact that the linear system ma- 
trix, which results applying the NEWTON’s method, is more simple than the 
corresponding one resulting from (3.4). 


e The CASH”s extended second derivative extended backward differentiation 
schemes [23] are particular cases of the following predictor-corrector schemes: 


k 
) ə OGUn+j— (ix, fanti- Rh? əə k> 
k+1 k+1 
) isinə k= m... bth? 27 Tifsi- k (3.6) 
One scheme step consists in the following stages: 


. compute ¥,, by (P); 
. compute 7,4, by (P); 


1 
2 
3. evaluate FGn41) and ITnai)3 
4. compute yn by (C); 

5 


. estimate error using the difference yn — Y,- 
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Moreover, the methods are chosen so that they be zero-stable and stable at 
infinity. The two classes studied by CASH are: 


k 
(P) Yn — Yn-1 = hh So Bibntink + hir, 
m i=0 
Class 1 (C) Yn —Yn-1 =h a + h?(yıga + Yeti Int): 
i=0 
k 
(P) 2593121: 
= 
Class 2 (C) Yn — axı =h a “haa. 
i=0 


e The SKELL-KONG”s formulae [53] consists in some combinations of the ADAMS’s 
type methods 


k 
(AM F+) “ln + 11 + h a — 0, 


i=1 


used in the integration of nonstiff problems, with the BDF, used in the integration 
of stiff problems, 


k 
(BDF) — ə QYnti-k thf, € 0, 
0 
in the form 


(AM FU") — sn (BDF) —0. 


3.1.4 Order 


The algebraic conditions of p order are the followings: 


k k k 
Sia? -q) Bi! q(q -1)) wit, = OS ap. 
go j=0 0 
More simplified, 


Pill) — joj) —- 20 — 196) =0, 2 ”1....p, 
Cp+1 = Ppt2(1) — (p + Lop(1) — plp + 196201) # 9, 
where 


k k 
201230: 51000: 


0 
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pan ($) = EE), oja lE) = Eole), öyə (6) = ESCE). 


The error constant is 
* Cpt 


o(1) 
k k 
= Tn Lot- (p+ 1) 97 A — (p+ DPA y 
j=0 0 


Lemma 8.1. 
(i) The ENRIGHT 5 formula with k steps has the order p = k + 2. 
(ii) The SDBDF with k steps has the order p ek +1. 
(iii) The SKELL-KONG’s formula with k steps has the order p=k +1 for any y™, 


(iv) If the predictor formula has the order p = k + 2 and the corrector formula, the 
order p = k+3, then the order of the CASH’s extended second derivative extended 
backward differentiation scheme is p= k +3. 


3.1.5 Stability 


The stability domain of a second derivative linear multistep formula is the set 


A={z€C| the polynomial p — zo — 276 satisfies the root condition}, 
where 5(€) = YE o yé. 


Im z 


-0.2-0.1 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 
Re z Re z 


Figure 3.1: Stability domains for the ENRIGHT’s formulae 


The characteristic equation is 


k 


X (ai — 28; — yi £” — 0, z — hA. 


0 
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For € = e, 0 € İ0,2zr) we get two root curves which describe the boundary of the 
stability domain. 

The DAHLQUIST’s second barrier is moved by two units. The greatest order of an 
A-stable second derivative linear multistep formula is four. This fact is a corollary of 
a more general theorem, which show us that the maximum order of an A-stable linear 
multistep formula is twofold than the maximum degree of the solution derivatives used 
in the formula (the DANTEL-MOORE s conjucture 1531). 

The ENRIGHT”s formulae break the DAHLQUIST’s barrier since the first two are 
A-stable with p > 2. The coefficients are given in Table 3.1, the stability angle in 
Table 3.2, and the stability domains in Figure 3.1. 


Table 3.1: ENRIGHT’s formula coefficients and the order of accuracy 


kl Pre Pros Pra Pr—s Pr—2 Prt Br Yk 


I 2 i j 
3 3 6 
i 5 29 -i | 
48 12 48 8 
zo -4 80 37 ç _ B] 
1080 20 40 540 180 
ır Lo sab df a, | 
5760 45 480 90 5760 32 
41 - 529 373 _ 1271 2837 317731 _ _ 863 | 
25200 40320 7560 10080 5040 604800 10080 
_ _ 731 179 _ 5771 8131 _ 13823 12079 247021 _ 275 | 
725760 20160 161280 90720 80640 20160 483840 3456 


Table 3.2: Stability angles of the ENRIGHT’s formulae 


kli 2 3 4 5 6 7 


o 190 90 87.88 82.03 73.10 51.05 37.61 
p13 45 6 7 8 9 


For solving the implicit equations we need some starting values. If we use an ex- 
plicit predictor formula, the scheme is equivalent to an explicit method (LAMBERT|74]). 
CASHİ25) has proved that the ENRIGHT’s onestep method can be applied in the fol- 


lowing scheme: 


(P) compute yo, yo) with the trapezoidal rule, starting with a value produced by 
the EULER’s explicit rule, 


(C) compute the approximate value of the solution at £, with the modified ENRIGHT’s 


formula 
52 yo 


6 


h 
Yn — Yn-1 = pu + fa-1) — 


where yn = yası — Yn, and the starting value yo, 
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This scheme is A-stable and of third order (an example for breaking the DAHLQUIST’s 
barrier). 
There are many improvements proposed for the ENRIGHT’s formulae: 


ə ENRIGHT [39] give up to an order unit in favour to solve more easy the linear sys- 
tem produced applying the simplified NEWTON’s iterations. The system matrix 
has the form 


(I — hök — hea, J”) = (I — hö,” > = —((,/2)” 


The other coefficients are determinated from the maximum order request. For 
k < 8 the condition system has two distinct solution, depending on Bp. The 
choise of one of these solutions is due to the numerical tests proposed by EN- 
RIGHT. 


e CHAKRAVARTI and KANEL [26] propose a modification in the idea to improve 
the order of the stiffly stable formulae. For stability at infinity, the ENRIGHT’s 
condition that o(x) = ypz" is successively replaced by the conditions 


a) o 
o) olz) = az (x + a)(a +b), 


The numerical tests on the influence of the parameters a, b, c prove the existence 
of some stiffly stable methods for 


(a) 3xkX8, (b) 3<k<10, (ce) 3<k< 12 


The two classes of CASH’s schemes [23] include some stiffly stable methods for 
k x 6. The maximum order for such a stiffly stable scheme is 9. The ENRIGHT’s 
methods are used only as predictor formulae. In Table 3.3 we have mention the 
stability angles. 


Table 3.3: Stability angles for the CASH’s scheme from the second class 


9 


23845 
45678 
o: 190 90 90 89 87 83 


The coefficients of the first four SDBDF are mentionated in Table 3.4. The stability 
regions are plotted in Figure 77. In Table 3.5 we mention the stability angles. 

The SKELL-KONG’s formulae have the same stability properties as the ENRIGHT’s 
formulae, but they are more easy to be implemented. The parameter y) is choosed 
to maximize the stability angle (see Table 3.6). 

The conclusions about the stability of the proposed particular cases of second 
derivative formulae are concentrated in the following lemma. 
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Lemma 3.2. 


(i) The ENRIGHT 5 formulae are A-stable for k < 2 and stiffly stable fork <7. The 
mazimum order of a stiffly stable formula is p = 8 (ENRIGHT /38/). 


(ii) The CASH’s schemes from the two class of extended second derivative extended 
formulae are A-stable fork <3 and stiffly stable fork <6. The maximum order 
of a stiffly stable formula is p= 9 (CASH [23]). 


(iii) The SDBDF are A-stable for k <3 and stiffly stable for k < 10. The maximum 
order of a stiffly stable formula is p= 11 (HAIRER, WANNER [53]/). 


(iv) The SKELL-KONG’s formulae are A-stable for k < 3 and stiffly stable fork < 11. 
The maximum order of a stiffly stable formula is p = 12 (SKELL-KONG [53]). 


Table 3.4: SDBDF coefficients for k < 4 


Ok-4 Qk-3 Qk-2 OAk-1 Ok Br Yk 


b | 


2 
1 
2 
1 
2 
1 
2 


0 —ü--.- 
—  49-10V147 
— 720 


-1 -0.8 -0.6 -0.4 -0.2 0 
Re(z) 


Figure 3.2: Stability domains for SDBDF 
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Table 3.5: Stability angles of the SDBDF 


[A] 1 2 4 5 6 7 8 9 10 


3 
2 3 4 5 6 7 8 9 10 11 
o: 190 90 90 89.36 86.35 80.82 72.53 60.71 43.39 12.34 


Table 3.6: Stability angles for the SKELL-KONG”s formulae 


10 11 12 


0 0.125 0.122 0.128 0.108 0.096 0.087 0.081 0.076 0.072 0.068 
90 90 90 89.42 86.97 82.94 77.43 70.22 60.68 47.63 28.68 


3.2 Multiderivative linear multistep formulae 


3.2.1 Formula 


A linear formula with k steps and m-derivatives is an iterative process of the form 


(3.7) 
where 2.70 020 2220 aşı, 77 0. Note l; = max{i| an 4 0). 
3.2.2 Particular cases 
e We get the standard linear multistep formulae when m = 1 and the second 
derivative formulae when m = 2. 
ə The onestep methods for k = 1 can be rewritten in the form 
i=1 ” 


For example, 


1. the ‘TAYLOR’s series correspond to a; = 0, 7=1,...,m, 


2. the formulae 


tn = Yaa +h S(—h) SO (yn) — PO (0) fO Cya) 
k=0 
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where m. 
P(x) = — L (a"(1 — a)" 
2... 
are A-stable (the PADE’s approximations to the exponential function 1591). 


e The ADAMS type methods are the formula with (aoi) = (1, 0, ..., 0, -1). The 
ENRIGHT’s formulae are an exemple of such methods. 


e The BRown’s formulae [12] generalize the backward differentiation formulae. 
The general formula is 


k m k m 
ə Oq ln4i-k = hi, LY, ə [a;l > 0, ə 19.) > 0, (3.9) 
i=0 j=l i=0 j=l 
with Bo = — Qk, (—1)+! 6; > 0, J — 0, — Dio oq = 0, Qk z 0, Br 7 0. 


3.2.3 Order 


A generalization of the barrier order is due to GENIN [48]: the maximum order of a 
k-step m-derivative formula is (m + 1)(& + 1) — 2. 
The BROWN’s methods have the maximum order p = k + m — 1 for 


—1)/ o. : 
25.5::5-55-5-- 
where i = 0,...,k — 1, 7 =0,...,m. 


3.2.4 Stability 


The order of an A-stable formula can be improved introducing the multiderivatives. 
In Table 3.7 we mention the order of such methods, function on the parameters k and 
m [68]. With * we have indicated a formula which is only A(a)-stable with a < 90°. 


Table 3.7: The order of the A-stable k-step m-derivative formulae 


The DAHLQUIST’s second barrier is generalized for the case of a multiderivative 
formula by the DANIEL-MOORE”s conjuncture. 


Theorem 3.1. (DANIEL-MOORE $ conjuncture [53]) The maximum order of an 
A-stable linear multistep formula with m derivatives is 2m. 
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The A-stability algebraic conditions have been studied by JELTSCH. Let 


k k (etl). 
— 9 uti, ri(z) = X ajz =(z-1) Pi z—l .—Ü,.... m, 
0 


0 
min(k,2) 1 m ok 
Qula Yo clear pO. Hod YOY aur". 
s=max(j,i-m) > i=0 7=0 


Proposition 3.1. 
(i) A multiderivative multistep formulae is A-stable iff (JELTSCH [63]) 


1. H(p,1) is an strictly HURWITZ polynomial, i.e it has all the roots in the set 
Rep <0; 
2. H(v,p) is a HURWITZ polynomial, i.e. it has all the roots in the set Rep < 
0, for anyv € R, 
(ii) A multiderivative multistep formulae is A-stable iff (JELTSCH /65/) 


1. r; are real nonzero HURWITZ polynomial; 
2. rii(p)/ri(p),? = 1,...,m, are real positive functions for p > 0; 


3. all the coefficients aj; have the same sign; 
(iti) A multiderivative multistep formulae is stiffly stable iff (JELTSCH [64]) 


1. it is Ag-stable; 


2. all the roots of the polynomial p(&)/(€ — 1) are inside the unitary complex 
circle; 


3. ifm; are the multiplicities of the roots of the polynomial p,(6) with |€| = 1, 
then miz, =m — j, 7 =0,..., miş 


4. if € is the root of pi(€) with (ŞI = 1, then all the roots of the polynomial 
Qmo(E, H) are positive real values, and if fi is the root of the polynomial 
Q mio($6, u) of k multiplicity, then 


Q. (6) “0 Vi, imici “ik mi, j=0,...,k-2. 
Application. Let 
2 1 2 1 2 
PoE) =E pilS)=—5(S +1), pal) = gU — 26 — 17), 


P(E) =- -E pE) =O + 26 + 1. 
Applying (iii) we can find that the formula is stiffly stable. 
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Lemma 3.3. 


(i) The BROWN’s formulae are Ao-stable and stiffly-stable for k < 3(m+1) (JELTSCH 


[62)), 


(ii) If the coefficients of the m derivative formula are independent from the stepsize 
h and ao; < 0, lal < Ong, q = 1,...,m, then the formula is A fm /(2m))-stable 


(Li [78]). 


3.3 Block formulae 


The building idea is to produce a set of approximations in different interval points 
at each integration step. The predictor-corrector schemes are the first examples of 
such composed methods. A RUNGE-KUTTA’s process can be also treated as a block 
methods, but we can evaluate the exact solution in some extra-division points. 

These formulae are also referred to as composed methods. 

The block methods have some strange properties: when we make a combination 
of some unstable methods we can get a stable block scheme or an greater order than 
those of the component formulae. These results was explained by the A-methods 
theory (ALBRECHT [2]). 

The block methods are suitable for parallel implementation. 


3.3.1 Formula 


The block methods lie in the iterative process 


(—1 i-1 
5 AijYknyi h ə Bi faq =9, 71=0,...,0-1, 


bn bn (3.10) 


where m is the number of starting points, / is the number of points at which we get 
the solution approximations at one step, and k is the advancing point number with 


kxi. 
Applying such a method, we follow the steps: 


[Step 0]: compute y1,...,Ym-1 by a starting procedure; 
[Step p]: follow the stages 


1, let £z; = Yrp4j, J = 0,... m — 1: solve the algebraic equations on the 
unknowns @m,---,;%m4l-13 


2. store the first output k values as the numerical values of the exact solution 
Y, Ykp+j = Uj, Where 7 = m,... m +k- l; 
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3.3.2 Particular cases 
e We get the standard linear multistep formulae for / = 1. 


ə The cyclic methods correspond to the case k = l. The set of / methods with m 
starting values is referred as a l-cyclic m-step method. The generic form as such 
formulae for the case k = m is due to BICKART and RUBIN [9]: 


AoY,, + A_1Yn-1 — h( BoY; + BuY..) — 0 (3.11) 
where Y, is the k-dimensional vector of all aproximations computed at one step, 
and Ao, Bo, Azı, B- are the [xl characteristic matrix with (A_1);; = (B-1)i; = 


0, Yj #0,7 =0,...,/—1. For each k, there is at least one k-cyclic k-step method 
with the order p = k [9]. 


e The VVILLIAMS- HOOQ”s formulae generalize the ADAMS methods: 
k 
Ynkr — Ynk+r-1 = h N Bri fakti» 1 < r < k. (3.12) 
j=0 
A more general formula is given in [123], referred to as advansed linear multistep 
formulae: 
k n+un 
N ajYnyj=k =h ə Bnj fn+j-k- 
j=0 J=Nn-Un 
e The BICKART-PICEL’s methods have the following form: 
l 
ə QijYkn+i T h fren; = 0, i= 0, ... „L. (3.13) 
j=0 
ə REINER [116] has studied some block implicit symmetrical formulae, defined on 


the symmetrical interval 1—1, 1] and a symmetrical interval division: 


l 
Yntj axı + Se Bin fn-14k: j3=0,...,8. 


k=0 
3.3.3 Similar methods 


Introducing the high derivative, we get the WATANABE implicit block multiderivative 
formulae [122]: 


Ynti = Yn th ə bb" fƏi, + bizh, 2(tn + 0:)h)), (3.14) 
j,k 
where z is the HERMTTE s interpolation function of the first p; — 1 derivatives of the 
exact solution y at t,4;. If Ọm are the HERMITE’s base functions, then: 


s pı—l 
2(ty + Oh) =X YO A Din (Oye 
İmzo m=0 


Note this formula with [po,...,ps]. The order of such a formula is at least 27020 pi — 1 
(WATANABE 11221). 
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3.3.4 Stability 
The DAHLQUIST’s first barrier is generalized by the following result. 


Theorem 3.2. (SLOATE, BICKART /121]) A block method which is zero-stable and 
stable at infinite has the order 


ps {mit 152, 


m 1, (—1. 
Moreover, if p> 1 andi = 1 then 


ps {™ k is odd , 


m+1, kis even. 


NEVANLINNA and SIPILA [91] prove the nonexistence of an A-stable block method 
in which all the approximations are explicitly computed. Moreover, the following 
theorem is a corollary of the DANIEL-MOORE’s conjuncture and a generalization of 
the DAHLQUIST’s second barrier. 


Theorem 3.3. /BICKART, JURY /7/) The order of an A-stable block formula can not 
exceed 21. 


Lemma 3.4. 


(i) The BICKART-PICEL 5 cyclic formulae are A(a)-stable for k < 10 (BICKART, 
PICEL [8/). 


(ii) The VVILLIAMS-HOOQ 5 formulae are A-stable for 1 < k < 8 (WILLIAMS/123/). 


(iti) The WATANABE’s following methods are A-stable: İp, pl and |p, p + 2] formulae 
for p> 1, and İp, p,p] formulae for p < 10 (WATANABE [122/). 


(iv) The REINER’s formulae have the order p > 2[(m + 2)/2] and are A-stable for 
1x5. 


For example, in Table 3.8 we mention the stability angle of the BICKART-PICEL 
formulae. 


Table 3.8: Stability angles for BICKART-PICEL’s formulae 


pek 


GEAR’s BDF 
BICKART-PICEL | 90. 
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3.4 Hybrid methods 


The linear multistep formula gives some approximate values of the exact solution 
at some points of the integration interval division. The RUNGE-KUTTA’s process 
produces some approximate values of the exact solution at some points between two 
consecutive points of the integration interval division. The hybrid methods is a combi- 
nation since it uses some evaluations at some division points and at some extradivision 
points. 


3.4.1 Formula 


A hybrid method with a single extradivision evaluation has the form 


k k 
33077 


J=0 J=0 (3.15) 


where frig = f (tn + Oh, yasə ), and yn4o is computed by a predictor formula. 


3.4.2 Particular cases 


e The ENGLAND’s methods [37] are implicit predictor-corrector schemes with a hy- 
brid corrector formula and with stability properties identically to the ones of the 
ENRIGHT’s second derivative formulae. The advantage of using these methods 
is the replacement of the second derivative evaluation with the approximation of 
the exact solution at an extradivision point. The general formula is the folowing: 


k 
(P) Haq = ə aj Yntj—k + ho f(ya), 


j=0 


k 
(C) Yn = Yn-1 + h S > Bi fat ik + hOf(Ynte)- 


j=0 


The predictor formula is referred to as the interpolation formula, and the cor- 
rector formula, as the quadrature formula. 


e The BEAUDET’s multidivision multistep formulae [5] are defined by the iterative 


process: 
m—-1 m-l n 
Ynti = X arp th SSS bri faitoo (3.16) 
0 j=0 kel 


where m is the stepnumber and n is the number of the used derivatives. 
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e The PATRICIO’s second derivative hybrid method [95] is the iterative process 


(P) Ynte + AYn = h(Boyl, + Boyse) + Poof + Yo frao) 


(C) yazı + doy, = h( Boy}, + Boule + AY a) +R oF, + Yo Fhe +7 fiy) 


One step lies in three stages: compute 7,,, by (P), evaluate yz, and Fao» then 
compute Yn4i by (C). 


e The CARROLL’s hybrid method [17] consists in the scheme 


On + A1Yn4o + O2Yn41 = Afati, Ynt = Yn + YAL — Ww) fn + fn), 
where 0< 80<1,0<u<l1. 


3.4.3 Stability 


The ENGLAND s hybrid methods have the same stability properties as the ENRIGHT”s 
second derivative formulae since the coefficients of the hybrid formulae are choosed 
such that the characteristic polynomials of both iterative process to be identically. 
Moreover, the accuracy order is the same. The request of k + 1 order for the interpo- 
lation formula and the request of k + 2 order for the quadrature formula implies that 
the predictor-corrector scheme has the order k + 2. For example, for k = 1 and p = 3 
we get the coefficients 


o(0) — 0(0 —1)) BO = —— 


ao(0) = (9 — 137, Bo(8) = ə o:(0) = —(9 — 2)0, 0:(0) = 30 = 2 


60 ” 6(9 —1) 
If 8 is choosed such that Bo = 0 then 
1 2 4 5 3 1 
Üzün qoy, qın =F Doz, Hr =F. 


The parameter 0 can be choosed also from the condition that the linear system matrix 
1 — h(Br + Bax) f’ — h?aß f? produced applying the NEWTON’s iterations for solving 
the implicit equations to be the square of some matrix (LAUTSCH İ761): 


2 
ağ = — (A) , 
2 
Lemma 3.5. 
(i) The ENGLAND 5 hybrid methods are A-stable for k < 2, stiffly stable fork < 7 
and of order p = k + 2 (ENGLAND /97)). 


(ii) The PATRICIO’s second derivative method is stable at infinity of Yoye = 
7072: Yo; Yı #0 and A-stable if 0 > Z, The A-stability coexists with the maxim 
order p= 5 (PATRICIO /95]). 


(iti) The CAROLL’s method of order p = 2 is L-stable (CARROLL /17İ). 
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3.5 Exponential fitting 


The basic idea of the exponential fitting method lies in building up a family of integra- 

tion formulae with some free parameters, others than the stepsize h. The restrictions 

imposed by the stability request do not affect the stepsize, only these free parameters. 
This method is very efficient in the transitory stage of the integration. 


Definition 3.1. A numerical method is exponentially fitted with the order s at a value 
hd if it exactly integrates, using a stepsize h, each initial value problem with a solution 
of the form P(t)e*’, where P is an arbitrary polynomial with the degree deg P < s. 


More precisely, if the approximation error in the component e” of the exact solution 
of the scalar test equation is 


T(z) = R(z)-e’, (3.17) 


and, for some value q = Ah, T(q) = 0 holds, then the numerical solution is exact in 
the discret meaning (yn = y(tn), Yn > 1 for the corresponding test equation) and the 
formula is exponentially fitted at q. The formula is exponentially fitted with the order 
s at q if q is a zero of the T function with s + 1 multiplicity. 

The exponential fitting may be accomplished by two techniques: 


e LINIGER-EHLE’s technique: determine the free parameters of a scheme in the 
purpose of the exponential fitting at some values (foe example, see CASH [21], 
JACKSON [61], LINIGER and WILLOUGHBY İSÜİ): 


e ISERLES’s technique: approximate the exact solution with a number of methods 
and compute the mean of the approximate values multiplied by some coefficients 
determinated by the exponential fitting conditions at some values (for example, 


see CASH [20]). 


The LINIGER-EHLE’s strategy for selecting the free parameters supposes that there 
is a group of eigenvalues of the JACOBIAN’s matrix near to the real axis. The formulae 
with a single free parameter are fitted to the mean of these eigenvalues. If there are 
two distinct groups of eigenvalues, we must use formulae with two free parameters. 
If the eigenvalues are changing during the integration, the formulae must be refitted. 
When the formula has many free parameters, one must be fitted in such a manner 
that the formula to be stable at infinity (exponential fitted at q = —oo). 

Example. Let the p-rule 


Yn = əz: = AU = u) fa + pfa]. 


For any value of u, the formula is exponentially fitted at zero with the order at least 
one. When the value u is choosed so that the formula is exponentially fitted to an 
arbitrary point q from the complex plane, then 


ula) = 1/q + 1/(e? — 1). (3.18) 


For u = 0 we get the EULER’s implicit rule which is exponentially fitted at zero and 
at —oo. The trapezoidal rule is exponentially fitted with second order at zero and we 
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get it for u = 1/2. The A-stability condition request that 0 < u < 1/2. This condition 
is not restrictive since u(q) : (—oo, 0] — [0,1/2] is an one-to-one map. The optimal u 
is the value for which 


— 1 — cd 
v = min max |R(q) — es) 


This optimal value is due to MIRANKER [86]: u = 0.122 for q = —8.19, and the 
objective function minimum value is v = 0.139. For EULER’s implicit rule the value 
0.204 is nearest to the optimal value than the value 1 of the trapezoidal rule. 


3.6 Extrapolation method 


3.6.1 Formula 


Let a sequence of pozitive integers ny < ng < ng <--- and the sequence hy > hg > 
h3 > +--+ related by h; = H/n;. Choosing the positive integers may be performed in 
various ways: 


1. the natural sequence: 1,2,3,4,5,6,7,...; 

2. the powers of the value two: 1,2,4,8,16,32,...; 

3. the powers 2” alternating with 1.5 x 2”: 1,2,3,4,6,8, 12, 16, 24,32,... 
4. the harmonic sequence: 2, 6, 10, 14, 22,34,50,.... 


Let a integration step-by-step method of p order. We evaluate, for each ? > 1, the 
numerical solution of the initial value problem on the interval with the length A, 
computing n; values at some points of an equidistant division with the stepsize h;. 
Note 

Tio = Yni- 


Let the interpolation polynomial 
P(h) = co + eph” + ei?” + ... + eppk-1 hP tT, P(h;) = To, 2 — J — k, ... 2 


We extrapolate the result by evaluating the limit value at h > 0: P(0) = eo = Tjk- 
The extrapolation table has the following form: 


Tho 
19 
: -Tipi (3.19) 
Tj-k+1,1 
T-ko 


The extrapolation method improves the initial approximations with each new co- 
lumn of the above table. 


Proposition 3.2. /53] The values T; p are the results of a numerical method of p + k 
order. 
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3.6.2 Stability 


The extrapolation of an A-stable method is not always an A-stable method. We 
mention the following counterexamples: 


the extrapolation of the trapezoidal rule with p = 1 and the natural sequence of 
positive integers leads to a non-A-stable method [53]. Moreover, when 


Tp) — Tiig 
j,k—1 j—1,k L k > 1, 
(nj/nj-k)}?— 1 


the resulting method of third order is not stable at infinity. 


Tjo =y; (H), Tix = Tic (3.20) 


CASH [19] proves that the extrapolation scheme based on the sequence {n;} = 
12, 3, 4, 6, 8, 12, 16, 24, 32, 48, ...), the AITKEN-NEVILLE’s iterative formulae 


T; k-1 — 1)-k-i L>] 
(nini) = 1 0 


and the EULER’s implicit rule, has an nonregular error behaviour in the stiff 


Tio = gÖtH)) Tik = Tik- + 


components. 


The A-stability conditions are satisfied when: 


1. we use some proper sequences nj; 

2. the approximations y,,(H) are processed before the extrapolation; 

3. we use some special basic methods. 

Examples: 

ə If we use only even numbers in the sequence nj, an interpolation polynomial 
with only even powers of the generating variable, the trapezoidal rule and the 
relationships (3.20), we get an A-stable method (HATRER, VVANNERİ53İ). 

e By smoothing the initial values we can get a method which is stable at infinity. 
Thus, let 

1 
Tho = qürəci + 2Yn, + Ynj+1) 
the trapezoidal rule and the natural sequence of positive intergers as stepnum- 
bers. The resulting method is L-stable (LINDBERG İ791). 
ə Some A-stable methods can be produced using the linear-implicit midlepoint 


rule. Let y’ = Ay + g(t, y) the differential system. The LAVVSON s function is 
y(t) = e“te(1). Transforming the initial system we get a new one in c: 


d= e~Atg(t, eAte), g(t,y) = f(t,y) — Ay. 


If we apply the midlepoint rule, 


—Afn Atn 


Cn), 


Cn41 = Cn—1 + 2he g(t, e 
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holds, thus, 
e yny = e” yn + 2hg(tn, Yn). 


We make the approximation et’4 œ~ ThA and we get the linear-implicit scheme 


(1 — hA)yı = yo + hg(to, yo), 


(I — hAjyg4ı = (I + hA)yr-1 + 2hg(to + kh, yg), k e1.....n) 


The method has the order p = 2. Suppose that we use the harmonic sequence. 
We make a smoothing step, 


İyo = (Yn, fə44 + Yn, /2-1)/2- 


and we use the iterative process (3.20). Then T;o, Ti1, T21 and Tz 2 are A-stable 
(BADER, DEUFLHARD [4]). 


The extrapolation scheme based on the EULER linear-implicit rule 


(I = hJ Masa — Yn) = hfn, Je Fy(tns Yn), (3.21) 


and the ATTKEN-NEVILLE s iterations leads to a method T;;, of the order k + 1. 
For the natural sequence we get A(a)-stability [35]. 


Chapter 4 


NONLINEAR AND 
MULTTVALUES METHODS 


The integration of some nonlinear systems gives us the idea to use some nonlinear 
interpolation formulae. The first question is concerning to the restrictions imposed 
in a stiff problem numerical integration. We can construct some explicit nonlinear 
formulae which are A-stable. 


Definition 4.1. A numerical step-by-step method is a linear one if applying it to the 
equation y’ = Ay, with A an arbitrary matrix, we get a linear equation in the discrete 
variable yp. 


The classical examples of linear methods are the linear multistep formulae, the 
predictor-corrector schemes, or the RUNGE-IKKUTTA’s process. 


4.1 Linear multistep formulae with matrix coeffi- 
cients 


The real difficulty in integrating a stiff problem is due to the fact that we must often 
compute some matrix inverses. 


4.1.1 Formula 


A linear multistep formula with matrix coeficients has the form (LAMBERT [75]): 


s k 
Or awa] 25 
i=l 0 


where Q, is a variable matrix choosed so that ||Q(n)|| is bounded and the matrix 


s—l 
22233 - farizeo (4.1) 
t=1 


J=0 


coefficient of y, is nonsingular. If bü) = Ü, 7 = 0,...,s — 1 the new approximation 
can be directly computed at each step with an unique matrix inversion. Practically, 
the —Q, matrix is an approximation of the JACOBIAN’s matrix. 
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The order of accuracy is independent on Qna. The method has order p if 


k mə 1 k 
(0 — m O — :m—17 (i) 
m =o day =u de! bi > 
J=0 0 0 
where m = .....p—, t=0,...,5. 


Some similar methods was proposed by LAMBERT[72] referred to as linear multistep 
formulae with medium variable coefficients: 


I 


k k 
Qj (tn )Yntj—k = h 2 Pia) fatit (4.2) 


where 
Aj(tn) =ajthaj(tr), Bilt) = Bi + hb;(ta) laO A, |bi(t)| < B. 
for example 


aj(t) = ajq(t), bid) = balt), a(t) = —(OF/Ay)(t, y(t)), 


MAJDA [82] (linear multistep formulae with dynamical filter), BJUREL [11] (modificed 
linear multistep formulae with exponential coefficients), NORSETT [93] (A-stable ge- 
neralized ADAMS-BASHFORTH method in which the coefficient of y,_1 is e””, where 


J = Of /Oy). 


4.1.2 Stability 


The advantage of these formulae is due to the fact that we can build-up some explicit 
A-stable formulae. For example, the second oder formula with s = 1 


(1450n) m —[ kora ağaz + (al +30) tna 


h 
= zl — @)fn-1 —(1+a)fn-2], —-l<ac<l. 
is A-stable. 


The behaviour of the classical linear multistep formulae at the scalar test equation 
is a prevision model of its behaviour when we integrate a linear system. This fact is 
not valid also for the formulae with matrix coefficients. 


Definition 4.2. A linear multistep formula with variable matrix coefficients is A- 
stable if all the solutions of the difference equation produced applying the method 
with Q, = —A to the problem y’ = Ay, with max, Re A; < 0, y(to) = yo, converge to 
zero when n — œ, for any stepsize h. 


Stability barrier: The maximum order of an A-stable formulae is 2s. 

Example. The method with s = 1, k=1, pl: a =1, ağ) = —İ, bit) =b, bO = 
1 — b, at) =q, ağ) = —a is A-stable if a +b > 1/2. For b = 0, we get an explicit 
formula, and for b = 1/2, a = 0, the trapezoidal rule. 
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4.2 Lambert’s nonlinear formulae 


Note that the EULER’s explicit rule is the result of the interpolation process with 
I(t) = At + B, Yn = T(t), Yn-1 = T(tn-1), fa — I'(ty-1) 


We get the iterative formula by the elimination of the parameters A and B. 


4.2.1 Formula 


LAMBERT [73] proposes the replacement of the polynomial / with a rational function. 
For example, we mention some conditions: 


dü) Mw = Mla) ton = Hai), fizi = Mba) 

ag o os ER mais Ihi) P= O12, fai = lor) 

a3 s FEF, gags = h) i= 01 fa = Pah Ha = P) 
The resulting formulae are the following ones: 

(L1) Yn — azi = Aade, 

gə) E 

(L3) Yn — Yn-1 = 7a 


These formulae are applied at each component of the solution. For example, for 


(L1) and the problem y’ = f(y), f =(f',...,f%), y(to) = yo 


: hal fi 1 
— ou = 1 1=1,..., N, möl 
Yn Yn-1 yta hfi ? ? ? Z 


The first two methods are examples of linear multistep methods with matrix coef- 
ficients: 


(L1) İT — h diag (Fa /Yn)] (Yn — Yn-1) = hfa; 


(L2) 1 + Lo Yn — hQnYn-1 + m + Lo Yn-1 = 2h fna—1, 


where —Q,, = diag 


2 Fii lyha ea] . 


t t 
Yn—1 7 Yn—2 
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4.2.2 Order and stability 


A disadvantage of these formulae is the locally unconsistency, i.e. the order of the 
method is not a constant (it is a variable function on the iterative values). This 
inconvenient can be avoid applying some linear multistep formula in a neighbourhood 
of an unconsistent point. 

The order p computed by LAMBERT [73] is 


> l, Yn- 20 > 2, fr—1 #9, > 2, fr #9 
PIL) £ 0, altfel, 74? — 0, altfel, PIS) = 1, fi = 0, fl £0. 


Lemma 4.1. /73/ 


(i) The formula (L1) is L-stable and its characteristic function is identically to the 
one of EULER $ implicit rule. 


(ii) The formulae (L2) and (L3) are A-stable and there characteristic functions are 
indentically to the one of the trapezoidal rule. 


4.3 Exponential methods 


Let the CAUCHY’s problem associated to the differential system y’ = Ay + g(t, y). 
A multistep exponential formulae is an iterative process of the form 


k k 
ə uu ə Öp, (Ah) Onti-k, ak 7£ Ü, (4.3) 


where a; are some constants independent on h, and öz,( Ah) are dependent on h and 
A. The formula is a nonlinear one, since the difference equation produced applying 
the method to a linear system is nonlinear in h. 

Practically, the exponentials are numerical approximated, for example by PADE’s 
rational function [53]. 

The advantage of such method is the possibility to use a reasonable stepsize, since 
the stability condition is ||®,,(AA)||Al|gy|| < 1. 

Examples: 


e The linear-implicit midlepoint rule defined in the extrapolation method section 
can be a first example. 


e For k = p= 1 the © function are (PREISER, LEE [77]) 


{ Piol Ah) = —(Ah)~?(a0(Ah — T)ei" — 1), 
®1,(Ah) = (Ah)? (aoet + I + Ah). 


e The generalized ADAMS-BASHFORTH method, due to CHU [28], has the form 


k 
aa = tty, +h ə Oi Yn—i41, Pki = 


i=0 
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1 k : 
—l 
-f e40-c)ə I] SFI- la, (4.4) 


0 Li J! 
For A = 0 we get the classical method. If the exponentials are approximated by 
the PADE’s function, the method is A-stable. 


e In the CERTAIN’s method one step consists in two stages 
(P) compute 


tn 
ye = efyn + f e74) gP (dt, (4.5) 
tn—1 


where g(t)? is the interpolation polynomial of k degree of the g(t, y(t)) 
values at tn-k-1,tn-k,---,tn-1. The output is the predictor value; 


(C) using the predictor value, compute 


tn 
5—:—. (4.6) 
tn—1 
starting from the interpolating points tn-k, tn-k+1;-- -tn 


If degg < k +1, then this explicit method is A-stable (MIRANKER [86]). 


4.4 Stiff separability 


In almost all cases, the eigenvalues of the JACOBIAN’s matrix of a nonlinear stiff 
differential system can be separated in two sets Mı, Mz, with the properties: 


(i) A, € Mi, Reà; € —1, |Imà;]| < əlReA,l: 
(ii) A, € Mə, |A;| İc, with p, c of moderate sizes. 


Definition 4.3. [10] The system x” = f(t, y) is stiffly separable, for a pair (to, yo), if 
the JACOBIAN’s matrix J = (Of/Oy)(to, yo) has the eigenvalues {A;} separable in two 
sets 

min |ReA;| > max İRe il 

1X:xs s+1<i<N 
The eigenvalues à1,..., Aş are referred to as the stiff eigenvalues. The eigenvectors 
associated with the stiff eigenvalues determine the dominant subspace, and the others, 


the dominate subspace. 


Let N the number of system equations and S = {1,2,...,N}. If the nonlinear 
system is stiffly separable, then there are two distinct subsets J and FE of 5 with the 
properties [UF = 5, IAN E = and 
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o(a) = {Aili € T}, (Jan) = {rele € E}, 
Re (AMY) < Re (AP?) <0, Wiel, e€ E 


The separability supposes that the subsets J and EF are not changing the dimension 
during the integration. This fact excludes a large set of stiff problems. 

The purpose of the separability is to reduce the computational effort by applying 
some explicit formulae to the nonstiff components and some implicit formulae to the 
stiff ones. 

Example. The HOFER’s method [57] consists in applying the modified midlepoint 
rule 


h _ 1 m h— 
Yn+1/2 = Yn + gin Unyi = Yn + h fa+1/2; Yn+1 = 5 (Ynti/2 T Unga F fmt): 


to the subset Æ and the implicit trapezoidal rule to the subset 7. The result is the 
following scheme: 


oa 


E h I I 
y = y) + stn ips = p + “Gö + foi): 


h — 
—(E E _(I I I Uu) 
(E) = yP + hfe igs (D = yi + TO + Fası/ə) 


E) 1, (Œ EŒ) , RE) D _ 
y = Ərə + mz + o Faai/a)) dü — nl 


An other example is the correction in the dominant subspace, due to ALFELD and 
LAMBERT [3]. 


4.5 (A,B)-methods 


The (A,B)-methods are also referred to as general linear methods, BUTCHER $ methods, 
or multivalues methods. 

These formula generalizes both the linear multistep formula and the RUNGE- 
KUTTA’s process. 


4.5.1 Formula 


A general linear formula is an iterative process of the following form 


(4.7) 


(7) 


i 


where y; ” are the external stages, and Y;, the internal stages. As usual, r = s inserting 


some zeros. Note u =r +s. 
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The matrix equation is 
y® = AyD) FAB fY), YO Ary") + hBif(Y), 


or 


YU) = AYD +4 hBF(Y™), 


where YU) = (y®, Y). 


We can rewrite the process in a symbolic form 


C Al Bı 
Aş bə 
An equivalent written form is that of an (A, B, C)-method. In the autonomous 
case, such a method can be described by the difference equation 


u u 


5x AHS 00100 as 
j=l 


j=l j=l 


The two forms are equivalents. The (A, B)-method is a particular case of a (A, B, C)- 
method, respectively when C = 0. If 


(n) ; 
—(n)_ J Yi i<u — {A 0 = [BC 
Yi = dz, >u’ us ə B= (4 0 ’ 
then the (A, B)-method produces the same approximations like the (A, B, C)-method. 


4.5.2 Particular cases 


e The (A,B)-method is explicit when B is an inferior triangular matrix, else the 
method is implicit. 


ə The RUNGE-KUTTA’s process with q stages is an (A, B) method with u = q+ 1. 
We must mention that there are some generalizations of the classical RUNGE- 
KUTTA’s process which are also examples of (A,B)-methods, like the multistep 
RUNGE-KUTTA methods, or the pseudo RUNGE-KUTTA methods. 


ə The explicit linear formula with k step, 


k k 
Yn = ə O kln-k +h ə Pk fn-k, 
j=l j=l 


can be write as an (A, B,C')-method with mal = Yn—i+1; U = k and 


ay wee OR] Ok Öl ... Dp 
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and as (A, B)-method with u = k + 1, yi, = Yn—k, and 


ay see Qk 0 0 Öl ... Dp 

1... ee 
45. po 

0Ü ... 1 0 0 0 ... 0 


If the formula is implicit, the first line of B must start with 60. 
e The predictor-corrector schemes can be also written as some (A,B)-methods. 


e Other examples of (A,B)-methods are the hybrid methods or the block methods. 


4.5.3 Order 


The order algebraic conditions was establish by BUTCHER [15]. 
Definition 4.4. The (A, B,C)-method is preconsistent if 


Ae=e (4.9) 
and consistent if it is preconsistent and, for some vector v (consistency vector) 
Av+(B4+C)e=v+te, (4.10) 
holds, where e = (1,..., 1)". 
Applications: 


ə a linear multistep formula is consistent iff 


k 


k k 
ıa-— od = ) Bi, v = (—1, —2,..., —k —1)7: 
j=l j=0 


j=l 


ə a predictor-corrector scheme is consistent iff the predictor and the corrector 
formulae are consistent; 


e acyclic scheme is consistent iff all its component formulae are consistent. 


4.5.4 Stability 


The stability of the (A, B)-method can be algebraic defined. 
Definition 4.5. 


(i) The matrix A is stable if there is a norm produced by a scalar product for which 


Aİ < 1. 
(ii) The (A, B)-method is stable if the matrix A is stable. 


Theorem 4.1. The (A, B)-method is stable iff the characteristic polynomial of the 
matrix A satisfies the root condition. 


4.5. (A,B)-METHODS 73 


The stability function is 
M(z) = Ag+ 2Bo(1 — zBi) A1. (4.11) 
A general linear method is 
(i) A-stable if p(M(z)) < 1, Vz € C_; 
(ii) Ao-stable if p(M(2)) <1, Vee R_; 
(Hi) stable at infinity if p( A — ByBy'A1) <1; 
(iv) zero-stable if p(A2) < 1. 


The DAHLQUIST’s second barrier is generalized by the following theorem. 
Theorem 4.2. /14] For an A-stable general linear method, r < s +1 holds. 


4.5.5 Convergence 
Let ||- || the euclidian norm on RY. The discrete values produced by a (A, B)-method 
are interpolated with an continous function y;(h), i =1,...,u. 


Definition 4.6. A general linear method is convergent if, for any autonomous 
CAUCHY’s problem with the standard condition for the unicity of the exact solution, 


2010012 a (=) wolf 0. ay 
holds for ? = 1,..., u. 
Let da real value, p a vector, and the norms 
N 
1 .0.5..54..2.x.. 


Definition 4.7. An (A,B)-method has p order of 
(i) convergence if Vk' > 0, 4K, R” so that: 


FO - YO in < k x YO -Yin < K, nh < Th <h 
(ii) consistency if there are w, k € N so that 
IFO - xO li < k (ATİ) -Yi < Ë, nh < T, Yh < h: 
(iii) stability if there is an a > 1 such that for any d and u > 1 


AN <a, [B"Alliya < game MBA < gae" 


Theorem 4.3. /31] 
(i) An (A,B)-method is convergent iff it is consistent and stable. 


(ii) An (A,B)-method has p order of convergence if it has p order of consistency and 
p order of stability. 
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4.6 Multiderivative multivalue method 


For the autonomous CAUCHY’s problem, the formula is [54] 


T s 


= yağlı” diyar af), 
k=1 


j=l j=l 


s 


1 n-1 k,1 — 
Y= Y ay E Y aG PED), . 
j=l kel jal (4.13) 


and the matrix formula 


ye) = Ayten + yy", HAOD (Y) 


(4.14) 


where A, A®,..., AC) are squared g + l-dimensional fixed matrix, characteristic for 
a method, q = r + s is the stagenumber, m is the maximum derivative degree, and 
ym — (y™® Y). 


Particular cases: 
e For m = 1 we get the (A,B)-methods. 


ə The explicit methods are characterized by some A inferior triangular matrix 
with ? = 1,..., m. 


ə The multiderivative multistep linear formulae are the particular cases when 


0 1 0 60 
008100 ... 0 "5 
A= . . : AC = : 
r 00.20 
1 pe? pe) 
Ü —d,ı —a I eto mil 


e The multiderivative RUNGE-KUTTA’s methods are some generalizations of the 
standard RUNGE-KUTTA process; in this case 


0... 01 a)... al) 0 
A= . . . AC) = . . : 
0 0 1 atı ak) 0 


4.7 A-method 


The class of A-methods include all the above mentionated formulae and schemes, both 
linear and nonlinear methods. 
The matrix formula is [2] 
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where Ö is the starting procedure, and A is a matrix which is independently on the 
solved problem. 
An A-method consist in 


1. a starting procedure, 
2. the advance formula, 


3. a correct validation 7(t,, h) of the approximate values Y;,: the global error is the 


difference Y, — Z (tn, h). 
Examples: 


ə The linear multistep formula 
k k 
2. oj Ynsi-k AX Bifa =0, =l 
j=0 0 


can be write as an A-method with Y, = (Yn,---,Yn-k), A= SOT, ® = (EON, 


where 
—Qk-1 —Qk—2 ... Tay — Qo 
1 0 Lee 0 0 
S= 0 1 0 0 ; (4.16) 
0 0 1 0 
k 2 
W(tn, Yah) = (1,0,...,0)"[Bef(tn +h, = XO YO + h(t, Yn h))+ 
j=l 


+X Bif(tn — (i 1)h, YY]. 


j=l 


e Anexample of modified linear multistep formula is the following difference equa- 


tion: 
Yn 1/2 1 1/2 1/8 Yn-1 
hy 1:1 0 0 0 0 hyi 
Yr1 | | 1/2 1/3 1/2 5/24 Yn-2 
hyni 0 l 0 0 hyn 


+h(3/8, l, —1/24, 0)7 O(tr-1, Yn-13 h), 


where @ is the solution of the equation 


1 1 
1: 
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If T'AT = diag {1, Az,.... Ar}, then [52] 


A=E | Az Fi, | ... A Er } E, 


where E = T diag {/,0,0,...}77', Fy = T diag {0, 7,0,..., Hİ, 


Let the local error 


enpi = Zlty + hsh) — AZ (ty, h) — h®(tn, Zlin, h), h). 
Definition 4.8. 


(i) An A-method has p order of consistency if, for a CAUCHY”s problem with an f 
function which is p times differentiable continuous, 


eo = O(h”), E(eo £ ei FE. + en) + enp = O(h?), 0 x nh € e, 


holds. 


(ii) An A-method is stable if A” is uniformely bounded for any n > 0. 


Theorem 4.4. /52/ Let an stable A-method for which 
eo = yo Fh + +++ + yp h? + OC), 
Ent = Yo(tn) + 1. ə: 


holds with y,(-) differentiable continuous. The method has p order of convergency iff 
has p order of consistency. 


Part II 


DEVELOPING NEW METHODS 
FOR STIFF PROBLEMS 


Chapter 5 
SPLIT FORMULAE 


5.1 Split ADAMS-MOULTON formulae (SAM) 


5.1.1 Motivation 


The classical ADAMS-MOULTON formulae, through implicit ones, can not be applied in 
the numerical integration of a stiff problem, since they stability domains are bounded. 
There are also the most used zero-stable methods for the nonstiff case. 

Specially for the integration of the stiff problems, CASH [24] proposes the split 
linear formulae. The split form of a linear multistep formula is 


k k 
ə OR-iYn—-i — h əə fa-i = hð f (ta Yn) + h( Be — 8) Fln, Yn) (5.1) 


where 6 is an arbitrary parameter, and the value yy is produced by a predictor formula. 
The formula class which was ”splitted” by CASH lies on the GEAR’s backward diffe- 
rentiation formulae. The result of this procedure is the improvement of the stability 
angle as for the GEAR”s formulae. In the numerical implementation, CASH have use 
one GEAR’s formula, and in the case of numerical unstability (due to some eigenvalues 
of the JACOBIAN’s matrix which are not inside the stability domain of the formula), 
it makes a switch to the corresponding procedure of the splitted formula. 
The study of the CASH”s idea generates two questions: 


Problem 1. Can we get some split linear multistep formulae for the integration of 
stiff systems, using some basic methods that are not efficient in the stiff case? 


Problem 2. Can we get an improvement of the stability angle also for other formulae 
than the GEAR’s ones? 


We will prove in the next section that the split ADAMS-MOULTON formulae pos- 
sess some stability properties which make them usefull for the integration of the stiff 
problems. Thus, it is possible to improve the CASH’s code. We can use an ADAMS- 
MOULTON’s formula as basic method and in the case of the numerical instability (that 
means we have detected an stiff problem) we switch it into the corresponding split 
formula. 
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5.1.2 Formula 


Note the ADAMS-MOULTON’s formula with the k steps 


k k 
Yn = Yai th SOF; VÝ fa yei dh). Bizi fais 


j=l j=l 


where 


and l l l 
Vf fj, VL = JT KT Sj- 


The local error of the k step formula is 
TERY = Fra h yE tn) + OCW), 


The characteristic polynomials are [99] 


pAM(6) = 76-10), gal) = PTE ED, 


wo- (E E E 
27 E 


t 


We propose the following predictor-corrector scheme with both implicit formulae: 


Predictor: Jn — Yn-1 = h ) (özi + BuO) fazi + (Or — BONE 


i=1 


k 
Corrector: ya — Yn-1 = h X Bp-ifn—i + (Br — O)hfn + ORF, 
i=l (5.2) 


where f,, = f (ta, y,), and Bp; are the coefficients of the classical ADAMS- MOULTON”s 
formula with k steps. 
In this form, the corrector formula is consistent and zero-stable. 


5.1.3 Order 


An efficient method in the implementation supposes that the both two nonlinear sys- 
tem of implicit equations, resulting from different equations, to have the same form: 
we get Op = —1. We use the simplified NEWTON’s iterations for solving the implicit 
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equations (an unique JACOBIAN’s matrix inversion at some start value). The matrix 
which is necessary to be inversed is the same for both difference equations: 


uu 0015:55—805) 
where yo is the starting value of the iterative process for solving the implicit predictor 
formula. 
The condition of a predictor formula of k order lies in 


VT9 = €], 
where v; = 2771 are the elements of some VANDERMONDE matrix, and 8 = 
(B,_13--+,89)?, ci = (1,0,...,0)7. Since V is inversable, the system has an unique 


solution B l l 
Bri = (—1 "Ch, 1=0,...,k. 
Thus, the split ADAMS-MOULTON formula with k steps can be described by the 


difference equations: 


k 
Predictor: Y, — Yn-1 = h ız. — (1Y C20) fni + (Be — OR Fn, 


i=1 


k 
Corrector: yn — Yn-1 =h X Br-ifn-i + (Br — Oh fa + öh). 
i=l (5.3) 


where (,_; are the coefficients of the ADAMS-MOULTON s formula with k steps. 
In [99] we propose the notation SAM for such a scheme. 


Proposition 5.1. The SAM scheme with k step has k +1 order. 
Proof. We use the identity 


k 


NO (=I Ckit = (= 1). 


0 
Since the classical ADAMS-MOULTON s formula has the order k + 1, the local error of 


the predictor formula is 


Tp — Oh ety FYE) 4 O(h**?), 


m 


Then the local error of the scheme is 


7 ðf R 
TEM = (Fpp yt (tn) un. nt? 4 One), 


where 7,41 is the error coefficient of the ADAMS-MOULTON’s formula. 
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5.1.4 Stability 
Characteristic equation 


The characteristic equation of the predictor formula is 


Psam(&) — q@sam(€) = 0, q=ha, 


where 


(5.4) 


For the corrector formula the characteristic equation has the two degree in q: 


psam(€) — qosam(€) + q”rsam($) = 0, q — hA 


The characteristic polynomials can be expresssed in terms of pam and Cam. 


Lemma 5.1. The characteristic polynomials of the split scheme (5.3) are 
psam(&) € Psam (£) = pam (E) = ETE- 1), 


osam(€) = (Öz — 20)pam(€) + o4m(), 
Tsam(€) = (Bx — 20)pam (E) + P(E- 1)". 
(5.5) 


Proof. Applying the scheme to the scalar test equation y’ = hd we get the difference 
equations 


k 
[1 — (br — *a\¥, = Yn-1 + 1) Oii — (-1)'C,O) yoni, 


k 
[1 = (Be — 9)alyn = Yn-1 +9 3.) Ba-iYn—i + Ons 


i=1 
where q = hà. We multiply the first equation with 1 — (3,8)q and we replace y, 
with its value from the first equation. It results an equation in second degree on q. 


Replacing y,_; > €”7” we get the characteristic equation. 


In terms of the (r,s) polynomials, the characteristic equation of the predictor 
formula is 


TSAMİZ) — @8sam(z) = 9, 


and that of the split scheme is 


rsam(Z) — qssamÜz) + @tsam(z) = 9, 
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where 


rsam(z) = Tsayv(z) = raniz) = ( + A , 


Ssam(z) = ram(z)— 9, 
SSAMÜZ) = (Br — 20)ram + sam (2), 
tsam(z) = sam(z) + 8°. 
Example. For the case k = 2, the characteristic equation is 


az? +bz+c=0, 


where 


2 5 1 
c = b(q, 0):— (1? + 50-5) q? (u 54 + 2. 


Studying the characteristic equations, we can prove the Ao-stability of some SAM 
schemes. 

Because the characteristic equation has two degree in q, the boundary of the sta- 
bility domain is defined by the following two curves: 


osam(e'") + y/o2 ale!) — 4psam(e)Tsam (ei) 


2TsAMİ€ ) 
osam(e'*) — voğa(e”) — 4psam (€ )rsam (€) 
2TsAMİ€ ) 


Ao-stability 


The proof that the split schemes are Ap-stable when k € 4 is not easy to be performed 
using the general formula. 


Case k = 1. The onestep corresponding method is based on the trapezoidal rule: 


(P) Yn — Yn-1 = Ohf (ta, Un) + (1 ~~ O)hf (tn-15 Yn—1)s 


(C) Yn — Yn-1 = Fina, Yn—1) + Ohf (tr, Yn) + G ~~ 0) hf (tn, Yn) (5.6) 


CASH [24] proves that this scheme is A-stable iff 0 > 1/4. 
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Case k = 2. The split scheme with two steps has the following form: 


(P) Yn — Yn- = a| (= -= 0) Sn Yn)+ 
+ E + 2 F(tn-15Yn—1) — (5 + 0) F(tn-2; 2] , 
(C) Yn — Yn-1 = 


5 2 1 
=h (= — 0) Fe. Un) + uu Yn-1) — pun. 2] + Ohf (ta: Yn) (5.7) 


where y, is the predicted value, and yn is the corrected value, the final one. 
The predictor formula has order two, and the split scheme, the order p = 3. If the 
system function is sufficient differentiable, then 


TE, = Oh? y) + O(h), 
and the local discretization error of the scheme is 
1 Of 


TE, = |—y — ?=(y, )y| ht hë). 
E Dy! Yə + O(h”) 


Lemma 5.2. /707/ If 9 € (—oo, —(1 + V6)/12) U [(v6 — 1)/12, 1/8], then the third 
order split ADAMS-MOULTON scheme (5.7), is Aov-stable. 


A strongest result was given in [99]: the (5.7) scheme is Ao-stable iff 


ZE “Hu 


. 5.8 
12 12 724 (5.8) 


For some -values the predictor formula has good stability property. For example, 
for 9 = —1/3, the predictor formula is A-stable, and for 0 < —1/12 is Ao-stable. 


Case k = 3. The split scheme with three steps is 


(P) Yn — Yn-1 = 1 (= = 0) Ftna Yn) + (= + v) F(tn—1, Yn-1)— 
? 30 t İ 0 t 
— (= + ) Ff (tn—25 Yn—2) + İs + ) f( — 
24 


(C) Yn — Yn-1 = 1 (= ~~ 0) F(tns Yn) + Silt Yn—1)— 


.. Un-2) + Z fh: 25) + Ohf (ta, Yn). (5.9) 


5.1. SPLIT ADAMS-MOULTON FORMULAE (SAM) 85 


The predictor formula has third order, and the split scheme, the order p = 4. The 
local discretization error of the predictor formula is 

TE, = 0hty + O(h”), 
and that of the split scheme 

19 


Of 
TE, = |—y®) — 9? 
720°" ay 


Lemma 5.3. /712/ If 9 € [-(2 + V10)/8, —(2 + V22)/24), then the split ADAMS- 
MOULTON scheme with three steps and of four order of accuracy, (5.9), is Ag-stable. 


(un )y | hë + O(hY). 


From the point of view of the stability and of the accuracy, the 6 optimal value is 
0 = —(2 + V22)/24. 


Case k = 4. Let the split scheme with four steps 
251 323 11 
P Tn — Yn-1 “hl 1-— —0 tn, y. — +46 tn—1,Un-1) — | — 
(P) ya — mi (= im Tn) + (+ im 1, Yn—1) (a 
+60 ) f(t )+ ag f(t )- P g f(t ) 
n—2: Yn-2 360 n—35 Yn—-3 720 n—4: Yn—4) |; 


(C) Yn — Yn-1 = 1 (= ~~ 0) F(tns Yn) + sitet Yn—1)— 


11 53 19 _ 
— ap t-> Yn—2) + 3601 -> Yn—3) — zən / (faza) ins) + Ohf (ta, Yn) (5.10) 


The predictor formula has four order, and the split scheme, the order p = 5. The local 
discretization error of the predictor formula is 


TE, = 0h>y®) + O(h), 
and that of the split scheme is 


3 ar 
— (6) — g2—— (5)İ pE h5. 
1607” mu | + O(h*) 


Lemma 5.4. /96/ For 0 = —(49 + 10V/147 )/720, the split ADAMS-MOULTON scheme 
of five order and four steps, (5.10), is Av-stable. 


Case k > 5. The split schemes are not Ao-stable. 


TE, = 


Proposition 5.2. For each k < 4 there are two values az, b, € R, so that the 
split ADAMS-MOULTON formula with k steps and k +1 order is Ag-stable at least for 
0E lak, dg]. 


Proof. Using the last lemmata, the request values are well defined: 


1/4, k=1 a, k=1 

— J —a, k=2 _ —(1 + V6)/12, k=2 
"k = $ (2 + /10)/8, k=3 "T ) -2+ V22)/24, k=3 
—(49 + 10\/147)/720, k=4 —(49 + 10V147)/720, k=4 


where a is an arbitrary real pozitive number. 
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Stiff stability 


It is posible that the proposed scheme are not only Ag-stable, but also stiffly stable. In 
this purpose, we study the form of stability domains. We have plotted the two curves 
which are defining the boundaries of the stability domains. 

In Figure 5.1 we have plotted the curves in a neighbourhood of the origin of the 
complex plane for the case k = 2 and some special values of the 6 parameter. 


Lop gc. 
[Ə = —(1 + V6)/12 — 


-0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.1 0 
Re(z) 


Figure 5.1: Root locus curves in the case k = 2 and for some 6 values 


Remarks: 


1. the boundary of the stability domain of the ADAMS-MOULTON”s formula with 
two steps corresponds to the case 6 = 0; 


2. for some 6 negative values, the stability domains includes the set C_, so the 
methods are not only stiffly stable, but also A-stable. For example, we have 
plotted the boundaries for 9 = —(1 + V6)/12 and —1/3. The methods are 
A-stable at least for 9 < —(1 + V6)/12. 


Because of the A-stability property in the condition of an third order of accuracy, 
the split scheme with two step is an example of breaking the DAHLQUIST’s second 
barrier. 


In the case k = 3 we get the plots from the Figure 5.2. 

For the value 0 = 0 the stability domain is inside the region delimited by the 
plotted curve (the classical method). For the negative 0 value, the stability domain 
includes the negative real axis; the method is not only Ag-stable but it is also stiffly 
stable. 

The constants from the stiff stability definition was computed for some 6 values 
using a numerical method (Table 5.1). 

The boundaries of stability domains of the split scheme in the case k = 4 and for 
two 6 values are plotted in Figure 5.3. For the negative 6 value the domain includes the 
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-1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 
Re(z) 


Figure 5.2: Root locus curves for split ADAMS-MOULTON scheme with & = 3 and for 


two 6 values 


real negative half-axis, i.e. the method is stiffly stable. The stiff stability parameters 
are D = 0.034, a = 88.260. The stability domain of the classical method corresponds 
to 6 —0. 


09 —n-..- 
— — 49-10V147 
~~ 720 
-1 -0.8 -0.6 -0.4 -0.2 0 
Re(z) 


Figure 5.3: Root locus curve for the split scheme with four steps 


5.2 Split second derivative multistep method 


5.2.1 Motivation 


The GEAR’s backward differentiation formulae are stiffly stable for k < 6. Splitting 
them, we get an improvement of the stability angles. On other hand, the GEAR’s 
formula with p = k = 3 is only stiffly stable, while the split scheme is A-stable. 
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Table 5.1: Stiff stability parameters 


—(2 + V10)/8 | 0.021 | 89.169 


CASH[24] suggests a numerical code based on the GEAR’s formulae and an algorithm to 
determine the numerical unstability and the switching procedure to the split schemes. 

The CASH’s idea can be extended to other formula classes than the backward 
differentiation methods. We mean some methods with stability properties requested 
by the integration of a stiff problem. For example, we consider 


e the ENRIGHT’s second derivative formulae (generalizations of the ADAMS’s for- 
mulae); 


e the second derivative backward differentiation formulae (generalizations of the 


GEAR’s formulae). 


5.2.2 Formula 


There are three natural ways to generalize the split linear multistep formula using the 
second derivative: 


1. splitting the first derivative, f, and the second derivative, g, of the solution 
computed in the current point of the division, tn; 


2. splitting only the second derivative of the solution, g, computed in the current 
point of the division, t,; 


3. splitting only the first derivative of the solution, f, computed in the current 
point of the division, t,. 


We analyse the results of the splitting method on two formula classes: 
e the second derivative multistep methods proposed by ENRIGHT; 
e the second derivative backward differentiation methods. 


We get some new schemes with stability properties requested by the stiff system only 
following the first two above mentioned ways. The stability angle can be improved 
only following the first way. 
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We propose the name of second derivative split linear multistep formula for the 
following generic equation 


k k k 
ə O k-iyn-i — h ə Bri f (tazi; yazi) +h? ə a-ig(tazı, yazı) = 


i=l i=1 


= N00 FCn, Ty) + h( 3x — 60) FCn Yn) — h20g (tn Ta) — R Ok — g(tas tn) (5.11) 


where y, is a value which is computed by a predictor formula, and g = y”. 
A particular case of the equation (5.11) corresponds to 6 = 0, and it is named split 
second derivative linear multistep formula: 


k 


k k 
ə Ak-iYn-i — h ə Bri fazi. ni) +h? ə Vig (tris Yn—i) + 


i=l 


+h?Og(tn, In) + h? (4x — 9) g(tns Yn) = 0. (5.12) 


with y„ produced by a predictor formula. 

For the convergence of the formulae (5.11) and (5.12), we use the coefficients 
Qk—is, Bk-i, Yk-i Of a consistent and zero-stable second derivative formula and a pre- 
dicted value produced by a consistent formula. 

We impose the following two rules for the predictor-corrector scheme: 


(i) we choose the predictor formula so that, for the case 9 = 0, the scheme is equivalent 
to the nonsplit corrector formula, i.e. the predictor formula has the following 
form 


k 
Onn — h(ök + BrO S nTn) +P? — OIl Tn) + > (ari + Thi) Yn i 


i=l 


k 


k 
—h ız. + By) fni Yni) +h? N Ori + 67; )G(tn—is Yn—i) = 0; (5.13) 


(ii) solving the implicit equations, we use a simplified NEWTON procedure with an 
unique matrix inversion and we suppose that the inversed matrix is the same for 
both predictor and corrector implicit difference equation, i.e. 


By = —6. 


If yo) is a starting value for the predictor formula, an initial approximation at 
t,, and we solve an autonomous problem, then the inversed matrix is 


29 A YO) = 7 h(6k — BO) Jy + h? (Or — 0)/7 


Of 
I — h( GB, — b0)— (y) +h 
(Br Emi. )+ Dy mill 


where J = OF /Oy(y). 


90 CHAPTER 5. SPLIT FORMULAE 
We resume us to the study of the order and stability properties of the schemes for 
which yoq =+- = yo = 0 and 
l. a, = —OK-1 = 1, Op-2 = ++» = ao = 0, By #0; 
2. ak #0, Öç = 1, Bk-1 = = Po = 0. 
Note that there are two cases: b #0 and b= 0. 


5.2.3 Formula with both split derivatives. Split ENRIGHT 
schemes (SES) 


Using as basic method an ENRIGHT’s second derivative formula or a second derivative 
backward differentiation formula, it is possible to improve the stability angle by the 
splitting procedure. 


Formula 


The proposed scheme, referred to as the split ENRIGHT scheme, has the following form: 


k 
(P) Tn = ei +h X (Bei + Fe) f Cni, Yni) + AlB — 69) FCn Ta) 


i=l 
k 


(C) Yn = Yn-1 + h > Br-if (tri, Un-i) + h( Bi, = bO) f (tn, Yn) 


+AbO F (tn, Tn) — R? (Ye — 9) g(t. Yn) — h?OgG(tns Tp): (5.14) 


where Bpi, yp are the ENRIGHT’s method coefficients, and 3,_; are computed from 
the maximum order conditions for the predictor formula. Note SES; the above scheme. 


Order 


The ENRIGHT’s classical formula with k steps is of k + 2 order. Note cz43, the error 
constant of this formula. 

Suppose that the predictor formula, in the given form, has the maximum order of 
accuracy, p = k + 1. From the order algebraic condition we get 


and the system 

VTB = €1, 
where v;; = i are the elements of a VANDERMONDE matrix, and 5 = 
(Br: 80), e: = (1,0,...,0)%. Since V is invertible, the system is uniquely 
solvable and 


By; = (10i fi,  i=1,...,k. (5.15) 
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By induction method, one can find that 


b= -54 (5.16) 


di 
Proposition 5.3. The split ENRIGHT scheme with k steps has k + 2 order. 


Proof. In the conditions (5.15) and (5.16), the predictor formula has the local trunca- 
tion error 


— SES 1 
y(tn) Un = TE, = man. + O(h*+3), 


so that, the SES scheme has the following local truncation error: 


d af 


TESES — R (k+3) th p —— 
n (ass (ta) + ha loy 


(ln: 0100 2550 


and its order is k + 2, the same like that of the classical ENRIGHT’s formula with k 
steps. 


Stability 


Since the ENRIGHT’s methods are A-stable for k = 1,2 (it is not necessary to improve 
the stability angle) our atention is concetrated on the cases when k > 3. Firstly, we 
study the possibility to obtain A-stable schemes, for k = 3. 


Case k = 3. The corresponding scheme is the following 


_ 307 11 _ 19 1 
25 th] (SE O) m1 + (= s+ 


3 7 1 19 2 — 
+50) F(tr—25 Yn—2) + İns — z0) 25 — (= + 0) h“ gltn, ya), 


307 11 19 1 
(C) Yn = Yn-1 + 1 (= + =!) Fe. Yn) + p-e Yn-1) — alu Yn—2)+ 
7 11 _ 19 : ; E 


The local truncation error has the following form 


3 110 
T B5PFs — (Su o eu) hÊ + O(h’). 


For the integration of a stiff system, we are asking stiff stability for the scheme 
(5.17). A necessary condition for this kind of stability is the Ao-stability. 


Proposition 5.4. /709/ For 0 > 0, the split ENRIGHT scheme (5.17) is Ao-stable. 
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-0.06-0.04-0.02 0 0.02 0.04 0.06 0 0.01 0.02 0,03 004 0.05 0.06 
) 


Re(hA) Re( 


Figure 5.4: Root locus curves for split ENRIGHT scheme with k = 3 and 220 


We have plotted in Figure 5.4 the root locus curves of the split scheme with k = 3. 
The curves are symmetrically to the real axis. To find stiff stability characteristics we 
analyse some domain surrounding the origin of the complex plane. For each specific 6, 
the stability domains include the half-planes defined by the plotted locus curves and 
the real negative half-axis. 

The plots show that, for some 0 values which satisfy the inequality 0 > 1/6, the 
split scheme is not only stiffly stable, but also A-stable and has a stability domain 
greatest than the one of the associated classical method. The A-stability property was 
tested also for smallest values than the above mentioned 6. For instance, for 0 = 1/12, 
the scheme is A-stable and has a resonable error constant which is dominated by ck+3. 

Note that the classical method corresponds to 0 = 0 and it is only stiffly stable. 
Thus the A-stable schemes break the DAHLQUIST’s barrier. 


Case k = 4. The corresponding split scheme consists in the following equations: 


3133 25 AT 41 
y = n— h — —é tn, Y, nn 40 tn- 9 AN — TOR 
3, =y al (Şa e) zu üzə m ++ 


7 
(P) +30) F(tn-2, Un-2) + (= = 50) F(tn-3, Yn—3) + ( m ət 
+70) f(tr—4y ins) — (5 + 0) h gftn, Ya), 


3133 25 47 41 
n — Yn- h Torn —0 tn, n — J (tn- »Yn-1) ə tn- s Yn— 
m = ma H| (SE +E) Hate) + Fma) -g ama) 


(C) + (a, Yna) — 


T 25 — 
45 F(tn-4; ins) = TIPE n 7n)— 


5760 
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32 


The local truncation error is 


13 
- (= R 0) h?gltn, yn) + Oh2G (tos Tp). (5.18) 


1 5 Of 
T EŞPE: — (Sau? Hə — 9” Yn 2 ht 4 O hè , 
864 12 Jy, ) A) 


so that the scheme is of order p = 6. Analogously to Proposition 5.4, we have proved 
that the above scheme is Ao-stable if 6 > 0. In addition, the scheme is stiffly stable 
for some values of 9. This conclusion derives from the analyse of the stability domain 
boundaries. For example, in Figure 5.5, there are plotted the root locus curves for 
some specific 9. 


-0.1 -0.075 -0.05 -0.025 0 — -0.170.160.15.0.140.13:0.120.11-0.1 
Re(22) Re(2A 


Figure 5.5: Root locus curves for the split ENRIGHT scheme with k = 4 steps and 
bz0 
From these curves we can see that the stability angle is improved when 6 increases. 


Note that the classical method corresponds to 0 = 0. In Table 5.2 we present the 
values of the stability angles for some @ values. 


Table 5.2: Stiff characteristics for SES with k = 4 and b 20 
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Remarks. 


1. If we pose the condition of the dominance of cz,43 in the coefficient error, then 
0 x 1/19. For those 0 small values, the improvement in the stability angle is 
not so large. For this reason, in the cases k > 5 the significant improvement of 
stability angle is possible only when the error constant is dominated by 6?. 


2. We expect that, at least for k < 7, the split ENRIGHT scheme with k steps to 
be stiffly stable for some positive real values 9 which depend on k. The proof of 
the Ao-stability of the general scheme for 6 > 0 is not easy to be performed and 
remains an open problem. 


5.2.4 Formula with both split derivatives. Split second 
derivative backward differentiation schemes (SSBDF) 


The proposed scheme referred to as the split second derivative backward differentiation 
scheme 


(P) OY + Yu (on i+ Opi) Yni 
+h( Bx — dO) f (tn, Gn) = hi — P)G(tns Ya), 


(C) OkUn T yuz ak if (tn is Yn öhu 
“(öh = OF (tn Yn) + BOF (nn) = Ke = OIC Yn) + KOIT) (5.19) 


where a-ı, Bk, Yk are the k-step second derivative backward differentiation formula 
coefficients, and @,_; are obtained from the maximum order conditions of the predictor 
formula. Note SSBDF;, the above scheme. 


Order 


The classical second derivative backward differentiation formula with k steps is of k+1 
order. Note with ck+2 the error constant of this formula. 

From the accuracy condition of the corrector formula, we get a value for the un- 
known variable 6. Supposing that the predictor formula has the maximum order k, 
the order algebraic conditions lie in the system 


Via = —2€2, 
where a) 
vij = m ” J Z 2, 
l, j=l, 
are the elements of an inversable matrix, @ = (@k-1;,---, Q0)", e2 = (0,1,0,...,0)7 
and 


-—.——. 


Since V is an inversable matrix, the system has an unique solution. The vector @ and 
the coefficient b are presented in Table 5.3 for each k which satisfies 3 < k < 6. 
The split scheme has the same order like the classical ENRIGHT’s formula. 
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Table 5.3: Coefficients @,_; of the predictor formula of k order 


_ _ 415 
822 822 137 822 822 

2806272 48208192 _ 23532768 24304896 _ 14624064 4862592 _ 3641472 

1524091 7620455 1524091 1524091 1524091 1524091 7620455 


Proposition 5.5. The split second derivative backward differentiation scheme with 
k steps has k + 1 order. 


Proof. The local truncation error of the predictor formula has the form 


SSBDF 


m 


TE Fon yO (t,) + O(h), 
Ak 


where ¢ can be expressed like a function on azı and 6 coefficients. The local trun- 
cation error of the split scheme is given by the following equation 


pt 
SSBDF _ 
TES = 


(crave + Ban wta) uu. 


Hence, the split scheme order is p = k + 1. 


Stability 


The classic second derivative backward differentiation method with k steps is A-stable 
for k = 1,2,3 and stiffly stable for k < 10. Hence, the improvement of the stability 
angle is important only in the case k > 4. 

Case k = 4. The split second derivative backward differentiation scheme with four 
steps has the following equations: 


415 92 3 159 4 84 
P) —-—7,4+ (44 791 — (24.79 1 t 124-789 | yn- 
(P) munun i HE h 1+ (F450) u 3 
1 17 25 42 1 
—~{—+—6) y,_ — — — 1 hfn Ya) =h? (= — 9) gliny 
m T 39 ) y “(ZF 30 ) F Un) (; ) 9g( Tn), 
415 3 4 1 25 42 
— n AYn—1 — ZYn- 99-37 1g$9n- “ə “mölh tn, n 
(C) rusu uru. (ÜR = 8) F(tns Yn)+ 


42 1 
"arena =I (2-90 at tin), (5.20) 
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The scheme order is five. The predictor order is only four: 


SSBDF. 144 83 
TE t = pp y(t, hê 
n IE To” (tn) + O(h’), 


2.00 
415 130 
The split scheme is Ag-stable at least for 9 < 0. The proof of this statement is similar 
with that of Proposition 5.5. 
We have plotted the root locus curves of the split scheme in Figure 5.6. The curves 
are symmetrically to the real axis. 


-0.05 -0.04 -0.03 -0.02 -0.01 0 
Re(22) 


Figure 5.6: Root locus curves for split second derivative backward differentiation 


scheme with k = 4 and b 40 


The curves show that, in the condition of Ap-stability, the 0 corresponding schemes 
are stiffly stable. Thus, the schemes with 6 < 0 are stiffly stable. In addition, for some 
0 values we can reach the A-stability property. We can note the increase of the stability 
domains with the decrease of the 6 values, and conclude that the scheme with 9 < —1/3 
are A-stable (breaking the DAHLQUIST’s barrier). 

If we put the condition that the error constant of the scheme to be dominated 
by the error constant of the classical method, then it results 191 < 0.245. For the 
maximum value, we get stiff-stability with an angle closed to 7/2. 


Case k = 5. The equations are 


120197, (5, 3799, 5, 7909, ; 
3600 ”” gə J tT 17 T xoş” py? 


(4 10 1031, 5 2491, (1,4, Q 
g 1ay 7 / "3 1167 gəz ” / 7 T 195 T xoş” j 
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60 137 2 
(C) ~~ — y + Əli ~~ Š pə + Hə ~~ a + sys + (=. 
3600 2 9 16 25 60 
225 225 2 


h 
— ——0 ) f (tr, yn) + OF (tn. 9,,) = (| — — 0 |) gftn, Yn) + h Oglta Yn), 
ER Hü Yn) + T37 F(a In) (5 al Yn) + h°Og(tns In) 


The scheme is of six order, with a predictor formula of five order: 


SSBDF. 3600 12019 
y(t, hi 
n 12019 7246607 Y CA TOA), 


and the local truncation error 


TE 


3600 11 12019 „ðf 
The scheme is Ao-stable for 9 < 0 and stiffly stable for some values of the 6 parameter. 
The root locus curves are plotted in Figure 5.7, to prove the statement of the stiff- 
stability. 


Figure 5.7: Root locus curves for the split second derivative backward differentiation 


formula with k = 5 and b 20 


The curves are symmetrically to the real axis. The stability domains are the 
regions defined by the root locus curves and the negativ real half-axis. One can see 
the improvement of the stability angle with the decrease of the 6 values. For 9 = —1/2, 
the stability domain of the split scheme corresponds to an A-stable method. 

The condition that the error constant be dominated by the error constant of the 
classical formula is not more restrictive than for the case k = 4: |0| < 0.221. There- 
fore, we suppose that the application of the splitting method to the second derivative 
backward differentiation formula leaves to an improvement of the stability angle not 
only for k = 4,5, but also for 6 < k < 10. 
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5.2.5 Formula with split second derivative 


By splitting only the second derivative from a second derivative backward differen- 
tiation formula, we get some stiffly stable methods, but not an improvement of the 
stability angle. 


Example 1. The split ENRIGHT scheme with three steps, 


307 19 5 1 
~~ IN h tn, Y 50 n=l; Yn— — an 
Un = Yn-1 + FEN wa) + (p Je 1: Yn—1) (+ 


uu + İns + z0) F(t n= sənə) m h? (= + 0) gü n Yn) 


307 19 
Yn = Yn- ün Fns Yn) + ggf -1 zi)” 


540 
-lra j+ üz: Hr G(tns Yn) +Oh7 gftn Jp) (5.21) 
20 n—2ə 1/n—2 1080 n—35 1/n—3 180 no Yn g ns Yn . . 


of five order, a predictor formula of third order, 


TE, = 0h*y® (t) + O(h), 
and the local error 


3 
800 


of 


y(t) +O? T 


X—ı yununu 
is Ao-stable at least for 9 € [—19/360, 19/360]. Plotting the root locus curves we get 


the information of the stiff stability for some 6 values (see Figure 5.8). 


-0.1 -0.075 -0.05 -0.025 0 -0.1 -0.075 -0.05 -0.025 0 
Re(22) Re(22) 


Figure 5.8: Root locus curves for split ENRIGHT scheme with k = 3 and b = 0 


The curves are symmetrically to the real axis. For each specific 6, the stability 
domains include the regions defined by the plotted locus curves and the real negative 
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half-axis (from the Ao-stability). Althought the proposed schemes are stiffly stable, 
the bigest stability domain corresponds to the case of the classical method, 9 = 0. 


Example 2. The split second derivative backward differentiation scheme with three 
steps, 


85 3 1 11 
P — —/ — 0).- — —20 141 — —0 ] yn-3 + — hf (tn, Yn) = 
(P) -n (B= Ona (4-20) (6-0 Şar 
-p (1-9 (tn, Yn) 
1 2 g ns Yn) 
85 3 1 11 
— pryn nm-17 7“ /n-— 019 —h İn, nj} = 
(C) xg Un + 3): — gaz + gUn- + hS Cn Yn) 
1 
of four order, with a predictor formula of second order 
TE, = 252102 + O(h*) 
and the local error 
36 | 1 Of 
TE, = — | yO (tr) — 20 (tn, y(tn) )yO(tn)| b” + ORS 
situ lu. 


is Ao-stable at least for 9 € 1—1, 1/21. 


0 0.02 0.04 0.06 0.06 0.1 
Re(22) 


Figure 5.9: Root locus curves for the split second derivative backward differentiation 


scheme with k = 3 and b= 0 


The plots of the root locus curves (Figure 5.9) for some 6 values show that the 
scheme is stiffly stable for 9 in the above mentioned interval. The classical formula 
corresponds to the case 0 = 0. 
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Chapter 6 


EXPONENTIAL FITTED AND 
HYBRID METHODS 


6.1 Parameter dependent formulae 


6.1.1 Motivation 


When we are applying a linear multistep method to a linear system of differential 
equations, y’ = Ay, the error en = y(tn) — Yn depends on the stability of the numerical 
amplification operator y,/Yn—1 and the closeness of this operator to the one of the exact 
solution. The numerical operator approximates the exact operator in a neighborhood 
of the origin. Note o(A) the spectrum of the matrix A. Then ho(A) is shrinking into 
a neighborhood of the origin by making A small enough. 

The idea of the exponential fitting method is to replaces the single point of the 
origin by a set of points in the complex plane, called the fitting points. The exponential 
fitting method becomes interesting for stiff systems when we observe that fitting points 
may be very large in magnitude so that h is not required to scale the entire spectrum 
of A into a neighborhood of the origin. 

We use the LINIGER-EHLE’s technique, i.e. derive integration formulae with free 
parameters, others than the stepsize, and then choose these parameters so that a given 
exponential function satisfies exactly the integration formula. 

The fitted formula are: 


e backward differentiation formulae; 
e second derivative backward differentiation formulae; 
e ENRIGHT’s second derivative formulae; 


ə LAMBERT’s nonlinear methods. 


101 
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6.1.2 Exponentially fitted backward differentiation formulae 
Formula 


We consider the following formulae depending on a real parameter a: 


k 


30 


10 (6.1) 


Order 
We get the maxim order p = k when a = Ü and the corresponding values a; are the 
coefficients of the GEAR’s formulae. In the case 
k 
X (a; + (=1)'Cha)yngj-k = hök fe. 
170 (6.2) 


where a;, Bp are coefficients of the GEAR’s formula and a Æ 0, the order of accuracy 
isk —1. Note BDF¢, the k step a dependent formula and BDF, the GEAR’s formula. 


Examples. We study the case k > 2, since the case k = 1 was studied by MAKULA 
et al. [83]: 


k=2 (3 + a)Yn — 2(2 + a)yn-1 + (1 + a)yn-2 = 2h fn, 


k=3 (11 +a)yn — 3(6 + a)yn-1 + 363 + a)yn-2 — (2 + a)yn-3 = 6h fn, 


k = 4 (25 +a)yn —4(12+a)yn-1 +6(6 +a)yn-2 — 4(4 + a)yn-3 + (3 + a)yn-4 = 12h fn, 


k=5 (137 + alın — 5(60 + a)yn-ı + 10(30 + a)yn-2 — 10(20 + a)yn-3+ 
+5(15 + a)yn-4 — (12 + a)yn-s = 60A fn, 


k=6 (147 + ayn — 6(60 + agacı + 15(30 + a) — 20(20 + a)yn-3+ 
+15(15 + a)yn-4 — 6(12 + a)yn-s5 + (10 + a)yn-6 = 60h fn. 
(6.3) 
Remark. For a = 0 we can get GEAR’s formula of order k, i.e. 
BDF? = BDF}. 


If a is chosen so that the last coefficient, depending on this value, is zero, we get the 
GEAR”s k — 1 step formula. Therefore 


BDF;! =BDF,, BDF3;?=BDF,, BDF; = BDF, 
BDF; = BDF, BDF;!° = BDFs. 


Thus, the proposed formulae are some generalization of the well-known GEAR’s me- 


thods. 
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Stability 


We know that the BDF are stiffly stable. We expect that the above formulae to have 
this property for other values of a as well. 


Proposition 6.1. /100] For 2 < k <6, the k-step backward differentiation formula 
of order k — 1, depending on a real parameter a is stiffly stable for any a € (cx, dx), 
where cp < di are two specific values. Moreover, for k = 2 and k = 3, there are two 
values c, and dpi, where ek < dk, so that the formula is A-stable if a € (ez, di). The 
constants are 


—2,k = 2, oo, k = 2, 

—5,k = 3, 13,k = 3, = = 
dy 4 -8,k=4, dp = k= = {ee y= {Or 

—24,k —5, 38,k = 5, .. O 

—15, k = 6, I,k =6, 


Im(z) 


0 —ü--.- 
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Figure 6.1: Boundaries of the stability domains for BDF and BDF 


Starting from the characteristic equation, we can get the curves which define the 
boundary of the stability domains. The curves are plotted for each k in Figure 6.1, 
6.2 and 6.3, and are symmetrically to the real axis. The stability domains lie in the 
exterior of the regions determinated by the plotted curves. 

Table 6.1 presents the stiff parameters for the above proposed formulae. 


Choosing the a value 


Studying the Table 6.0 we can remark that the optimal stability domain is obtained 
at each k for the smallest a negative value. On an other hand, the error constant of 
the k step k—1 order formula depends on the a value. For some small absolute values, 
the error constant is smallest than the corresponding constant of the classical formula 
with k — 1 steps and k — 1 order. 

Choosing the optimal value for the free parameter, in both sense of stability and 
accuracy, it is necessary a compromise between the two tendencies: an |a| large value 
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Figure 6.2: Boundaries of the stability domains for BDF? and BDF? 
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Figure 6.3: Boundaries of the stability domains for BDF% 


from the stability condition and a small lal value from the order condition. Our 
solutions are mentionated in Table 6.2. 

The free parameter a can be used for exponential fitting at some exponential func- 
tion. The fitted methods are also stiffly-stable. 


Proposition 6.2. /708/ For any2 < k <5, the backward differentiation formula with 
k steps and k—1 order, exponentially fitted to an arbitrary negative real value, is stiffly 
stable. Moreover, for k = 2 and k = 3 the formulae are A-stable. 
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Table 6.1: Stiff stability parameters for BDF% 


If the test equation is y’ = Ag, A € R7, and the used stepsize is h, then 


2q — 3)e"4 + 4e? —1 
ə -i nü 
(6q — 11)c”7 + 18c” — 9e4 +2 k 
q 3 ? — 3, 

a(q) = 7 (6.4) 
(12q — 25)e*4 + 48e? — 36e” + 16e — 3 L4 
000 ey o”? m 
(60q — 137)e€?7 + 300637 — 300e*! + 200c?7 — 75e? + 12 has 

(e! — 1) ? 
where q = AA. 
Remarks: 


1. The k-step exponentially fitted formula at q = 0 is GEAR”s formula of order k. 


2. The exponential fitting to q = —oo reduces the k-step formula to one with k —1 
steps. 

3. We assume that the six-step exponentially fitted formula is also stiffly stable. 
The numerical calculations are very difficult and tedious, so we omitted this 


proof. 


106 CHAPTER 6. EXPONENTIAL FITTED AND HYBRID METHODS 


6.1.3 Exponentially fitted second derivative backward diffe- 
rentiation formulae 


Formula 


We propose the following iterative process as an a dependent second derivative back- 
ward differentiation formula: 


k 
X (a; + AO; )Un+j—k + h( pr + abp) fr = hye On 


170 (6.5) 
where g = df /dt. 


Order 


We get the maximum order p = k + 1 for a = 0 and a;, Öp, the coefficients of the 
second derivative backward differentiation formula with k steps. In the case 


X (a; + (D(a — DİL“ (k — 1) /? yma je + hl Be — a) fe = hegn: 


3-0 (6.6) 


where oy, Bk, yr are the coefficients of the classical formula, the order of accuracy is 
k. Note SBDF¢, the a-dependent k-step second derivative backward differentiation 
formula and SBDF;, the classical second derivative backward differentiation formula. 


Examples. Let the following a dependent formulae: 


6a — 7 2a —1 3 — 2a 1 
k=2 Yn + 2(1 — a)Yn—1 + —— Yn-2 + hfn = =h gn, 
4 2 2 
66a — 85 3(2a — 1) 1 —3a 11 — 6a 1 
k=3 —— ou, ə(1— n— — ~ Yn- ——1o- hfa = =h? gn 
gg W FBC = a)yn-1 + —— Yn- + — Yn- + Gh fn = hg 
paa 300a = 45 a aya p D, p A380 
— 144 Yn @)Yn-1 9 Yn-2 9 Yn-3 
4a —1 25 — 12a 1 
——— Yn- ——h. n — nh n 
+ in Yn—4 + D f g 


Remark. For each k and a = 0 we can get the classical formula with k steps: 
SBDF? = SBDF;. 


If a is chosen so that the last coefficient, depending on this value, is zero, we get the 
classical formula with k — 1 steps and k order: 


SBDFİ” = SBDF), SBDFİ” = SBDF,, SBDF!* = SBDF3. 
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Stability 


The SBDF;, formulae are stiffly stable for k < 10. We can try to find some conditions 
on the parameter a which conserve this property. 


Proposition 6.3. /700/ For2<k <4, the k-step second derivative backward diffe- 
rentiation formula of k order, depending on a real parameter a is stiffly stable for any 
a € ap, where 


1 k=2, 
a, = l 14/15, k — 8, 
5/6, k —4. 


Remark. We assume that, for k < 10, any k-order k-step second derivative backward 
differentiation formula is stiffly stable for some a values. 


The boundaries of the stability domains are plotted in Figures 6.4 and 6.5. The 
curves are symmetrically to the real axis. 
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Figure 6.4: Boundaries of the stability domains for SBDF? and SBDFS 


The stiff stability parameters are mentionated in Table 6.3. 


Table 6.3: Stiff stability parameters for SBDF% 


0 
1/8 
1/4 
1/2 
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Figure 6.5: Boundaries of the stability domains for SBDF4 


Choosing the a values 


Similar to the case of the backward differentiation formulae, we can get an improvement 
of the stability parameters when a increases and the local error descreases when |a| 
descreases. Thus, it is justified to use the & step formula in the place of the classical 
formula with & — 1 steps, if a is chosen like as in Table 6.4. 


Table 6.4: The optimal a values for SBDF% 


k | Order | Error(SBDF%) . our “a -stability 
pek — Co , di) 


ol 2 (-1/2,1/2) 
31 3 (-1/3,1/3) (0, Mts) 
Al 4 (-25/12,25/12) -00, [1/8.5/6) 


The existence of a free parameter in the condition of stiff stability make possible 
to apply the exponential fitting method. We expect that the property of stiff stability 
is extended to the fitted formulae. 


Proposition 6.4. /108] /f2 < k <4, then the second derivative backward differentia- 
tion formula with k steps and k order, exponentially fitted to an arbitrary real negative 
value, is stiffly stable. 
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If the test equation is y’ = Ag, A € R7, and the used stepsize is h, then 


1 (24? — 3q + 7)e"4 — Set +1 k=? 
2 (3 — 2q )ey —4er1 ? .. 
_ J 1(189” — 66q + 85)e*! — 108e74 + 72e! — 4 L 

ala) = 1677 Tfl — 6q)e” 8e 498-2 =3, (68) 
1 (7247 — 300q + 415)e*! — 576e?! + 216e”! — 64e! +9 L4 
12 (25 — 12q)e™ — 48e”! + 36e” — 126? +3 no 

where q = AA. 
Remarks: 


1. The k-step exponentially fitted formula at q = 0 is the second derivative back- 
ward differentiation formula of order p= k + 1. 


2. The exponential fitting to q = —oo reduces the k-step formula to one with k —1 
steps and with the same order. 


3. We assume that, for any 3 < k < 10, the exponentially fitted formulae are also 
stiffly stable. 


6.1.4 Exponentially fitted ENRIGHT second derivative formu- 
lae 


Formula 


We propose the a dependent ENRIGHT second derivative formula: 


k 


Yn = Yn-1 th X (B; + 4B, ng ik +? 16g 
J= (6.9) 


where g = df /dt. 


Order 


We get the maximum order p = k + 2 for a = 0 and 0), yk, the coefficient of the 
ENRIGHT’s formula. In the case 


k 
1 
Yn nı + h S (6; = (—1)' Cha) Yui + haaa. 
0 


(6.10) 


where 6,, yk are the coeficcients of the ENRIGHT’s formula, the order is k + 1. Note 
ESDF¢ the k step a dependent formula, and ESDF;, the ENRIGHT”s formula with k 
steps. 
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Examples: 


2 1 1 
k=1 m=nath| (Z-a) f+ (540) faa] = (2-2) a 
29 3 5 1 1 
k=2 Yn = Yn-1 +h (2 = Za) fa + (= +20) fazi = (z + xe) fes —- 
1 
~~ (t-a) Nh? dn, 


307 11 19 103 
k= n = Yn- h Tn o n TR n-1— | 57 re n— 
0 -.: (= re) fot (ta) l (z+) at 


(6.11) 


Remark. For a = 0 we can get the ENRIGHT’s formula of order £ + 2: 
ESDF?. = ESDF;,. 


If a is chosen so that the last coefficient depending on this value is zero, we get the 
ENRIGHT”s formula with k — 1 steps and order k — 1: 


ESDF,'/"*=ESDF,, ESDFŞ”/””” = ESDF). 


Stability 


The ENRIGHT’s formulae are stiffly stable for k < 7. The question is whether the 
formula with variable coefficients is also stiffly stable. 


Proposition 6.5. /100] For 1 < k < 3, the k +1 order k-step ENRIGHT second 
derivative formula, depending on a real parameter a, is stiffly stable for any a € ag, 
where 


1/6 , kİ, 
a, = € 1/24, k =2, 
1/144, k =3. 


The boundaries of the stability domains are plotted in Figures 6.6 and 6.7. The 
stability domains include the unbounded sets of complex points from the left half-plane 
limited by the plotted curves. 

The stiff stability parameters are dependent on the a parameter (see Table 6.5). 

Remark. We assume that the formula with k-steps and a free parameter a is stiffly 
stable also if 4 € k x 7. 
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Figure 6.6: Boundaries of the stability domains for ESDF? and ESDFS 
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Figure 6.7: Boundaries of the stability domains for ESDF% 


Choosing an a value 


The stiff stability parameters are optimal for the smallest negative a values and the 
local error decreases with the decrease of the |a| value. A compromise can be proposed 
(see Table 6.6). 

The free parameter a can be use by the exponential fitting method. We expect 
that the exponentially fitted formulae are also stiffly stable. 


Proposition 6.6. /708/ For 1 < k < 3, the ENRIGHT second derivative formula of 
k+1 order and with k steps, exponentially fitted to an arbitrary negative real value, is 
stiffly stable. 
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Table 6.5: Stiff stability parameters for ESDFZ 


a} P |e fale} ele Do) 


1/12 | 0.086 | 86.241 | 1/48 | 0.077] 87.850] 1/288 | 0.162 | 87.010 
0 | 0 0.098 | 87.922 
-1/6] 0 
-1/3] 0 
-1 | 0 


fe Error(ESDF¢) <Error(ESDF;_1) | Stiff stability 


—00, dı) 


2 (-1/24,1/24) (-00,0] 
3 (-7/360,7 /360) (-00,-7/720] 


If the test equation is y’ = Ag, A € R7, and the used stepsize is h, then 


(q? — 4q + 6)e! — 2q — 6 


k=1 
6lq(q— Let +q] 
_ J (6q7 — 29¢ + 48)e74 — (20q + 48)e* +q B 
e(a) = 24[q(2q — 3)e” + 4qc1 — q] ? 52 (6.12) 
(11447 — 614q + 1080)€”7 — (513q + 1080)e”? — 54qe? — Tq 
og 7 owen 6 04 kg 


180[q(6q — 11)e” + 18qe*? — 9qef + 2q] 


where q = AA. 
Remarks: 


1. The k-step exponentially fitted formula at q = 0 is the ENRIGHT’s formula of 
order k + 2. 


2. The exponential fitting to q = —oo reduces the k-step formula to one with k —1 
steps. 


6.1.5 Exponentially fitted multistep nonlinear formulae 


Suppose that the initial value problem function has an autonomous form. 
The LAMBERT’s nonlinear formulae are built up using some rational interpolation 
functions, instead of the classical polynomial interpolation function. 
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The error problems encoutered by the nonlinear methods can be diminuated using 
the exponential fitting method. We mention some results in this direction. 


Example 1. Let the rational function 
At+B 
(3 ——”. 
t+C 
and the conditions for an iterative process which modelates the solution of the initial 
value problem 


Yn = T(tn), Yn-1 = T(tn—1), fa — T (ta—1) 


If we eliminate the coefficients t, B and C, we get an equation in the unknows 
Yn-1; Yn, Jn-1 and A. Supposing that the numerical solution is exactly integrating 
the equation y’ = Ay, i.e. Yn = e”^yn—1, we get the value A and the iterative process 


(6.13) 
The method is an onestep one and is of order 


> l, Yn- 0 
nE I 17 


= 0, otherwise. 
For q + —oo the f,-1 coefficient is one and we get the first LAMBERT’s formula, (L1). 
A most important property is decribed in the following proposition. 


Proposition 6.7. [103] The formula (6.13) is strongly A-stable for any q € R”. 


Note that the formula is an explicit one and it is A-stable, i.e. it breaks the 
DAHLQUIST’s second barrier. 


Example 2. Let the rational function 


At+ B 


I(t) = —” 
(2) E+Ct+D 


and the conditions 
Yn-i — T(tn-i), ? — Ü, 1,2, fn-1 — Tu). 


Eliminating t, B, C, D we get an equation on xə, Yn—-1, Yn, fn-1 and A. Let 
A = chy,_1. Supposing that the numerical solution is exactly for a scalar test equation, 
we get the iterative process 


(Yn—1 — Yn-2 (ch fazi — Yn—1) + 2hayy-ə fn-1 


Yn = Yn-1 F Yn-1 > ————, 
i i YZ + 2€Yn—1(Yn—1 — Yn—2) — h.fa-i (eyazı + 2Yn-2) (6.14) 


where 
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The method is a two step one and of the order 


= 0, Yn-1 = fa =0 
P2) > ə altfel. 
For extremely large q values c œ% 0 since limj4—. c(q) = 0. 
Proposition 6.8. /96/ The formula (6.14) is strongly A-stable for any q € R”. 
Example 3. Let the conditions 
At+B 
I(t) = = 
0) P4Ct+D 
Yn—i = T(ty-i), ? = 0, l, fn-1 = I'(ty-1), f'(tn-1) = P (ty-1). 


Supplimentary, we suppose that the numerical solution is exact for a test scalar equa- 
tion. Then 


Yn = Yn—-1t 
+hy (2fn-1 + hff_1)¥n-1 + 2h(e — 1) 44 
n=-13 3 Ty POOP DOP 
2yn—1 t+ hlel fazi - hfa) = (faz +h fia) na + 2h? fai (6.15) 
where Det L 9fet—1 
c  c(q) :- ag = 2) + 2e = 1) 
q(q — 2)e4 + q(q + 2) 


This is an onestep method and has the order of accuracy 


= IL Fa — 0, 
PS) > 9 altfel. 
Two special cases can be distinguished. For the first one, we get by exponential fitting 
at —oo (stability at infinity), since lim, c(q) = 0. The second one is a LAMBERT’s 
formula, (£3): we get it by exponential fitting at zero. 


The exponential fitting to an arbitrary negative real value q lies also to a strongly 
A-stable formula. 


Example 4. Let the rational function 


AVP + Bt+C 
j= ———€” 
#?4+ Dt+ek 


the conditions 
Yn-i = I(tn—i), i = 0,1,2, fna-5 = L(ta-;), j=1,2 


and the request of the identity between the numerical solution and the exact solution 
of the equation q = hà. Then 


Un = Yn-1(Yn-1 — Yn-2)(4cyn-1 — (€ — lJyn-2) + 2h(1 — eye fa-2— 
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—hYyn—2(Yn—2 + CYn—1) fn—1 = 2h? cyni fn-1 fn-2, 
Un = (3c + 1)Yn—1(Yn—1 — Yn—2) + 2h(1 — C)Yn—1fin—2 — h(Yn—2 + CYn—1) fa” 
—2h? f.—ı fn-2, 
eq — (24? —q 11) (24 E1)e qa 
(q — 3)e° + (2q + 7) (242 +3q Ele 


The formula has maximum third order. The exponential fitting at —oo lies to the 
formula 


where 


e74 
C = 
e24 


(Yn—1 = Un-2)) + h(Yn-1 ~~ Un-2)fa-ı = 2h? fn-1fn-2 


n = un 6.17 
” 3 : 2(Yn—1 nü Yn—2) ~~ h(2fr—2 nü fn—1) ( ) 
Example 5. A formula with good implementation results is the following: 
_ (Yn—1 — Yn—2) Fr 
Yn = Yn-1 + 2......... “yr,” 
al — (Yn-1 — Yn-2) fra (6.18) 


The formula is a two step one and has maximum third order. It is exactly for any test 
equation y’ = Ag, if the first iteration produced by the starting procedure is exactly: 
note Wn = Yn/Yn—-1, we get 


(1 = TARLA 


1 “ii. 
7 i M — (1— wae)? 


= Wn-1 


From this equality we get the A-stability property for the proposed formula. 


6.2 Onestep hybrid methods 


6.2.1 Motivation 


Introducing some extradivision solution evaluations in a linear multistep formula, the 
hybrid methods allow the breaking of DAHLQUIST’s second barrier. The onestep hybrid 
methods are closed to the form of RUNGE-NUTTA”s methods. The actual variants do 
not request a large calculus volume, but, for high stability properties, are supposing 
that the extradivision values depend on yn. 

A particular class of such methods is due to ENGLAND[37]. The basic idea is to 
maintain the same stability properties like for ENRIGHT’s second derivative formulae: 
the second derivative evaluation is replaced by an evaluation of the solution in an 
extradivision point. 

Note that, for the onestep case, the ENGLAND’s formula of third order has the 
following form: 


39 —1 39 — 2 h 
mu — n h—— n h—— n — Barn 1391 ? 
Ynti = Yn + i —og fa + 60 = 1)" +1 6000 — 1) +6 


2-25: 1 
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One can get the particular formula studied by ENGLAND from the condition of a zero 
coefficient for f, in the first equation: 


h 3h 
Maxi = Yn t+ qün + “q Haki fa: 
(6.19) 
5 4 Ph 
Yn+1/3 = g yer + gin — dott 
The local discretization error is: 
ht Of 
TE, = —— ly (t,) —4— (tes yta yP tn Oh). 
ən U (tn) ass... 


The scheme is L-stable. 

The hybrid methods are some particular cases of pseudo RUNGE-KUTTA methods 
(for details see [96]). For example, the BOKHOVEN’s methods [119] are represented by 
the iterative process: 


Yn+1 = yn +hY_ biki, ki = (1 = 9i)un + liyny +h X aiki, i= 1.....s. 


i=l j=l 


In this 6-class, we distinguish one A-stable scheme of four order: 


h 2h 
nd = Yn + ian + fa) + FEL 


il h (6.20) 
Ynti = Yasi + Yn) — g — Fa). 
The local discretization error is: 
h5 of 
TE, — — OG) —5— (tn, y(tn) yO (tn O(h® 
go V Ch A Yta) yi (tn)| + OCR?) 
Another A-stable scheme, of third order, is the folowing: 
h 
Yn+1 = Yn T o faq + İns), 
2-3 4 — v3 h 
Un, = 6 Un T 6 Ynt1 — glint (6.21) 
4 — v3 24+ /3 h 
Unto = 6 Un T 6 Unt + gin 


where 9; = (3 — V3)/6, 6, = (3 + V3)/6, and the local discretization error 


ht Ə? 
2101107 


We mention also the OBRECHKOFF’s A-stable formula of four order, with second 


derivative 
h h? , , 
Yn+1 = Ün + fe + fa) = To (aza = fi), 
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and the local discretization error: 
h 


TE, = — 
720” 


(tn) + ORS). 


In the following sections we study three formula classes, generalizing the above 
mentionated scheme in the idea to get some A-stable or L-stable scheme of superior 
order, with a similar calculus effort. 


6.2.2 The first method class 


Formula 


A generalizing scheme for the ENGLAND’s 6-class and the BOCKHOVEN’s method 
(6.20) is the folowing: 


v+u — 
+1 -—i J Jangi T fn) — hv fna40, 
O0—w++r 0—w-r 
40 = WYn 1 — w)yn + h——— fn h— — fn 
6 = WYn41 + (1 — w)y z Jai + z S (6.22) 
which has the minimum order in each equation. 
We discuss also the following variant: 
l+u l-—wu 
Ynt1 = Yn + — ni (vYn+ı + (1 — v)yn) + ah fata, 
h h (6.23) 
Una = WYn41 1 (1 w)Yn f -00 ulu £) frat + 5 (8 — tü — zf. 


2 


The first equation is referred to as the quadrature formula, and the second, as the 
interpolation formula. 

We search the order and the stability properties to find some new methods which 
have atmost the same performances like the methods mentionated in the last section. 


Order 


Since, in the given form, the order is p = 1, the scheme is consistent. 

The maximum order for the scheme (6.22) is p = 4. The corresponding scheme is 
the BOKHOVEN’s method of four order which has the same stability function like that 
of the OBRECHKOFF’s formula. 

The maximum order for the scheme (6.23) is p = 3: we request the third order for 
the quadrature formula and the minimum second order for the interpolation formula. 
We get two formula classes with a free parameter in the interpolation formula. The 
fixed coefficients are, for the first class, 
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respectively, for the second class, 


1 2 1 


where w is the free parameter. If the interpolation formula has only second order, the 
local discretization error of this formula has influence on the local discretization error 
of the scheme. If we request the third order for the interpolation formula, the local 
discretization error has the same error coefficient like the quadrature formula. From 
the two formula corresponding to the condition of third order for the interpolation 
formula, only one is A-stable. For this the 6 value is 0 = Z and the scheme has the 
following form: 


h 3h 
Yn+1 = Yn + qin + yala 
əəə A Aq 
Yn4+2/3 = əy hai oz ün D7 n+1 D7 ne 


The local discretization error is 


(6.24) 


4 
~ 216" 


Note that the error constant, 1/216, is lowest than the ones of the ENRIGHT’s second 
derivative formula, and of the ENGLAND’s scheme. 


TE, Dn) + OCh). 


Stability 


Let the case of a third order quadrature formula and a second order interpolation 
formula. The conditions of strong A-stability restrict the w parameter: 


- 2/3, v=0 
€? 1/3, v=1' 


For example, the scheme (6.24) is strongly A-stable and has a third order interpolation 
formula. 
Exponential fitting and stability at infinity 


The free parameter w can be used in the purpose of the exponential fitting. If the 
fitting point is q € R_, then 


4 (q? — 6)ef +2 +q4+6 


2M77774174170 — 0 
(q) 9 (q? — 2q)le +g H2 7 07 
w = 
1 1 (q? + 6q — 24)et + 547 + 18q + 24 : 
— — ——--—..---..---..-...--..---...------- v= il. 
9 (q? — 2q)€1 +q? + 2q ? 


The schemes are strongly A-stable for any q € R... We get some particular schemes 
by exponential fitting at zero and infinity. The corresponding coefficients are the 
followings: 


lim w(q) = 


q——oo 


(ln v=0 lim w(q) = 


2/3, 0-0 
5/9, v=’ q—ü 


1/3, 0-17 
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The exponential fitting at —oo is equivalent with the request of stability at infinity. 
We get two formulae. The first one, for w = 5/9, is the ENGLAND’s onestep method. 
The second one, for 0 = 2/3, w = 8/9, is the folowing: 


h 3h 
Yn+1 = Ün + Ir + ru. 
8 1 öh (6.25) 
Yn+2/3 = goH + gin — dnt: 


The local discretization error is 
he 


TE, = -— 
216 |” 


OYE.) + ıı y(tn))y(tn)| + OR’). 


The error constant is comparable with those of the ENRIGHT’s second derivative for- 
mula and the corresponding ENGLAND’s onestep method. The scheme has also the 
L-stability property. 


6.2.3 The second method class 
Formula 


The BOKHOVEN’s scheme uses two additional points. The following method, with 
two extradivision evaluations and of the minimum order one, is a generalization of the 
BOCKHOVEN s scheme: 


+u 1—xu 


fni +h fantz 


F (L = wi)yn 


F (1— w2)yn 


Order 


The maximum order of a such scheme is four. The corresponding method is the 
following: 


h h 
azi = Yn + ZI n+1/24V3/6 + J n+1/2-V3/6 

fl 2v3 1: 2V3 hf1 V3 
 .xnuı.. LS Fo pee GL EG Init (6.27) 


hf1 V3 
tg (=$) Jn 


The local discretization error is 


hə 5 əs , i 
TEn = T30 G (tn) + Lot e) + O(h®). 
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Comparing the error constant with those of the BOKHOVEN’s method, and the 
OBRECHKOFF’s formula, we ascertain a substantial improvement. 

If we request only a third order for the scheme, we get three free parameters. Let 
u, W 1, wg the free parameters. Hence, 


1 VƏ — üz?) 1 30 — uz) 1--2u 2u—1 
bi = = 1, b = = = —_ _, 25 = — 7°. = — ., 
2 6(1 +u) 2 6(1 — u) 6(1 +u) 6(1 — u) 
Stability 


The stability function of the method of minimum third order is 


14+ (1 — öt)z + (1/3 — t)z? l+u l-—u 
= — — — — t= 
He) 1 — 21: — (1/6 — ():2 ” q Ty 


w2 


The condition of strong A-stability requests that £ > 1/4. The BOKHOVEN’s third 
order scheme is only A-stable, since £ = 1/4. 

If t = 1/3, we get a method subclass which depends on two parameters (w1, w2). 
All schemes from this subclass have the same stability function, therefore the same 
stability properties like the ENRIGHT’s onestep second derivative formula. 

If w, = w = 2/3, then u is the only one free parameter which can be chosen so 
that the local discretization error has a minimum constant error. 

The stability function of the four order scheme is identically with that of the 
OBRECHKOFF’s formula, and, therefore, the scheme is A-stable. 

In the class of third order schemes we can distinguish some L-stable methods. If 
we request a third order in both equation of the scheme, then one solution has the 
property of the stability at infinity: 


2+ V3 2— v3 
Ynt1 = Yn + 4 —— hf 4.V33 + hn 5/3" 


9 — 2/3 2/3 3 — v3 wie 6 (6.28) 


Yny 3/3 = əə + Yn = Bh Int + Sh 
9423 2/3 34+ V3 Wee 346 
Yn—/3/3 = ~g Vu mü m in = 9 hfnti — — əş YT. 


The local discretization error is: 


4 
TE, = Py (ty) + O(h*). 
72 

Note the identity of the constant error of the above scheme and that of the ENRIGHT’s 
onestep second derivative formula. The derivative evaluation is replaced by two func- 
tion evaluations. The advantage of this scheme compared with the BOKHOVEN’s third 
order scheme consists in the nonexistance of a perturbation in the local discretization 
formula due to the interpolation formula. 
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6.2.4 The third method class 
Formula 


The BOKHOVEN’s third order scheme uses three function evaluations at one step of 
the NEWTON’s procedure for solving the implicit equations in the unknowns yn4e,, 
Yn+0.> Ynti- We can estimate the number of function evaluations to 3m + 1, where 
m is the medium number of the NEWTON’s iterations to one step of the integration 
interval. Usually, the value m is small, since the convergence order of the NEWTON’s 
iterations is quadratic. 

Let the following scheme which also needs three function evaluations at one NEW- 
TON’s iteration: 


+1 = Yn + ahf(uyngr + (1 — uyn) + öh F(Ynta)+ 
tech f(vyngi + (1 — olma), 


+0 = Ulan + (1 m W)Yn + hdf (uyn4a + (1 m DUR 
--hef(oyazı + (1 — olm). (6.29) 


Order 


The conditions of third order are the folowings: 


1 
a+b+c=1, @=wt+d+te, 0? = wt 2(du + ev), au + 00 + cv = 5, 


1 1 
au’ + DY + ev" = 5, au + BO" + cv = 5. 


Therefore, the third order scheme class depends on three free parameters. Let u, v, 6 
these free parameters. 

The local discretization error is the sum of the error produced by the quadrature 
formula, 


hi 1/1 302 — 40 +1 Pf 
2.80 


and the error produced by the interpolation formula (in the approximation of the exact 


value y(t, + 9h)), 


TR) = F o — by) — (epee i.a ur (24) n) ım 


“Tale. y(tn)) + Oh’). 


The ENGLAND’s onestep method is a particular case of the above mentionated 
schemes, which has the third order. The maximum order for a class scheme is four, 
corresponding to the BOKHOVEN’s onestep method. 
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Stability 


Getting some L-stable methods, firstly we request the property of stability at infinity 
for our scheme. Algebraically, this condition means that 


d(l—wu)+e(1—v)=0 
Hence, the method coefficients are 


y—t vy — $ ung u— 0 1 


ou OP 0)0v -u)” ~~ une! e u o) T RY 


v—u u — v 


where 6 is the solution of the equation 


[6uv — 3(u + v) 4319? — [6uv — 2(u + v) + 2]0 + w = 0. 


Note that we get a subclass which depends on two parameters, u and v. All the sub- 
class schemes are L-stable, since the stability function is identically with the stability 
function of the ENRIGHT’s second derivative onestep formula. 


Examples. 
ə Foru =0, v= 1l or u = 1, v = 0, ve get the ENGLAND s 6-class. 


ə The values v = 1, 0 = i lie to the ENGLAND’s onestep particular method. 


correspond to a method subclass which depends on the 


wiw 


ə The values v = 0, 0 = 
parameter u: 


3h h 
Untl = Yn + TP+ + I 
1 2h 2(1 — u)h 


8 
Yn+2/3 = goH + oi” — ulu + (1— uyn) 4 


For u = 1 we get the scheme (6.25). 


ə Ifw+v=1, then the coefficients expressions can be simplified: 


1 1 
pla ge ty 
2 6 6(u — v) 6(u — v) 6 

9—1 9—1 
a= 2 b=1, c= —. 
v — 4 u—v 


The local discretization error of a such formula is 


TE, = -2 | + və (ou) + 6u(1 — u) (r) (tn, 0) - 
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E (; - £) gəbə (FEF) iste + OU) 


For simpliffing this expressions, we can request that 1 — 6u(1 — u) = 0. The 
resulting scheme has the form 


6 


W =m HAS (va) FEAS 
hf (G — Pyn + G +) 


(4 + o) Yn— (6.31) 


For some system functions, the constant error of this formula can be lowest 
than the one of the ENRIGHT’s onestep formula, which has the same order and 
stability properties. 
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Chapter 7 
PARALLEL COMPUTING 


7.1 Motivation 


The problem associated with the stiff ordinary differential equation systems in parallel 
processing is due to the fact that the calculus can not be started simultaneously on 
many processors with an explicit formula. We propose two different algorithms. One 
is based on the method of extrapolation which has stability properties requested by 
stiff equations. This algorithm can be applied to any kind of stiff systems. The second 
one is built-up for some special classes of stiff ODE, y/(t) = A(t)y(t) + g(t). The both 
algorithms have a high efficiency when the system function has many components. 
The approximation error is equal to the one produced by the analogous sequential 
algorithm. 

Parallelism in solving ODE can be expreseed via three distinct avenues (ISERLES, 
NORSETT [60]): 


1. coding a specific method so that it can be performed simultaneously on several 
processors; 


2. splitting variables in a multivariable ODE system between processors; 


3. exploiting parallelism in performing the request computer algebra, solving linear 
and nonlinear algebraic systems of equations. 


7.1.1 First avenue 


Many serial methods contain a potential degree of parallelism. Most suitable for the 
parallel implementation are the multivalue methods. The parallel methods from this 
avenue preserve the stability and accuracy of the basic sequential algorithm. 

The predictor-corrector schemes can be easily implemented in a parallel mode. 
Many exemples are presented by MIRANEER and LINIGER ([87], [88]). The application 
of those schemes to stiff system is not succesfully because almost all ones are equivalent 
with some explicit formulas. The stability characteristics can be improved. In this 
direction we can mention the papers of KATZ et al. [70], WORLAND [125], GHOSHAL 
et al. [49]. These improvements make possible to handle slightly stiff equations. 


125 
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The block methods are easily adapted to a parallel mode with no degradation in 
the accuracy of the solution. The performance of HUTCHINSON and KHALAF parallel 
implementation [58] is dependent on the number of scheme nodes. CHU and HAMIL- 
TON’s algorithm [29] has a numerical efficiency depending on the dimension of the 
solving system and on the complexity of the system function. 

ISERLES and NORSETT [60] investigate the degree of the parallelism of the RUNGE- 
KUTTA methods. The parallelism depends on the exploitation of the sparsity structure 
of the RUNGE-KUTTA matrix. KARAKASHIAN and RUST’s algorithm [69] is designed 
to solve a linear system of ODE by a RUNGE-KUTTA process. The efficiency of this al- 
gorithm depends on the system dimension. For the proposed test system with variable 
dimension, the parallel mode becames competitive only for a number of components 
of hundreds order. 

The parallel RUNGE-KUTTA method proposed by EVANS and SANUGI [42] makes 
use of the special form of the method matrix. In an other paper, EVANS and MEGSON 
[41] consider a systolic array construction of extrapolation tables. The basis schemes 
are not suitable for the stiff integration. 

GALLIGANI and RUGGIERO [44] propose a parallel method which is suitable to be 
applied to a large set of linear ODE systems. The method has good stability properties, 
which makes it usefull for the stiff systems integration. 


7.1.2 Second avenue 


This way request a reconstruction of the class of numerical methods for ODE. The 
basic idea is due to NIVERGELT [92]. The integration interval is divided into equal 
subintervals and each processor is responsable for the integration on only one subinter- 
val. To make efficiently the algorithm, the processors must be started at appropriate 
time and must work simultaneously. For this purpose it is used a starting value (an 
approximation of the solution at the begining of the processor subinterval). In the 
special case of a stiff system, the starting value can not be obtained by an explicit 
method, only in the case when the subintervals number and the step size are very 
small. In the case of using an implicit starting formula, we reach the class of block 
methods. An other problem associated with the stiff case is due to the fact that the 
approximate solution, produced by the parallel algorithm, can be unstable, also when 
the basic method for integration on each subinterval has good stability properties. The 
error of the approximate parallel solution is different from that of sequential solution. 


7.1.3 Third avenue 


KNIRCH [71] studies a scheme for solving the stages algebraic equations of a RUNGE- 
KUTTA method. The reason for unsuccesfully using of this scheme in the stiff case is 
the same like for NIVELGELT method. 
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7.1.4 Proposed algorithms 


In the frame of the first avenue, one algorithm explotates the parallelism in the extra- 
polation methods. It can be applied to any stiff system for which the basic extrapo- 
lation method works. The efficiency depends on the dimension of the system, while 
the degree of the parallelism is function of the requested accuracy of the numerical 
solution. 

The second algorithm is an hybrid one. It can be applied to linear stiff system with 
variable coefficients. The proposed parallel algorithm divides the responsabilities of 
the processors on the interval of integration (second avenue), produces the same error 
like the sequential algorithm (characteristic for the first avenue) and makes use of the 
linearity of the system (third avenue). The efficiency depends on the complexity of 
the system function, ans the degree of parallelism, on the dimension of the system. 

In both cases, the numerical results show that a high efficiency can be obtained for 
a reasonable system dimension. 


7.2 Basic concepts on the algorithm efficiency 


Knowing an algorithm for parallel computing, we note: 
T, - the execution time of an implementation of the algorithm using p processors; 


Tı - the execution time of the same implementation of the algorithm using only one 
processor; 


To - the execution time of the best implementation of the most closed sequential 
algorithm. 


For the given parallel algoritms the efficency may be computed in many ways: 


Enum - the numerical efficency of the parallel algorithm is computed by the following 
formula 


Lo 


Enum = nn." 
Tı 


(7.1) 


ET - the efficency of the parallel implementation of the algorithm is given by the 
next realationship 


Ep - the efficiency of the parallel algoritm with p processors is given by the formula 


The ideal values of these efficiency measurements are 1. The practical values are lowest 
because: 
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ə for Enum, almost all sequential algorithms are very difficult to divide in a number 
of equal units (in the sense of the same computing effort); 


e for He", the communication time between the processors is greater that the 
computing time necessary for an arithmetic operation. The ratio is between 500 
and 1000, depending on the network; 


e for Ep, since Enum < 1 and Bl <1. 


The notion of the speedup, Sp, is related to the second measurement: 


Naturally, the maximum value of 5, is p. 
The employing factor of a specified processor, U, is the fraction of the execution 
time when the processor is busy. 


7.3 First algorithm 


7.3.1 Problem to be solved 


We consider a stiff system in the general form 


yi) = F(t y(t), € 0/7), 


with the initial condition y(0) = y°. Note with D the dimension of this system. The 
integration interval is divided in N equal subintervals of H lenght. 


7.3.2 Basic method 


We consider the extrapolation method proposed by DEUFLHARD [35], based on the 
linearly implicit Euler rule: 


(I= hJ (yiti — yi) = hf lti yə), 17> 0, 


where J is an approximation to (Of /Oy)(to, y(to)) and yo = y(to). We take a sequence 
of integer numbers, representing some step numbers, ny € mə € ---, respectively by 
n; =1, i > 1. We define 

Tio = yn, (to + H), 


the numerical solution obtained by performing n; steps with size h; = H/n;. The 
values are extrapolated according to 


Tik-1 — 1)-1k-i 


Ti = Tip 
gk jk t+ n/n; —1 


kl, 


so that Tj, represents a value given by a method of k + 1 order. To the next step, we 
take 
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and to > to + H. 
DEUFLHARD [35] proves that T;; is A(a)-stable, at least for 7 < 8. Thus, the 
method has stability properties requested by the stiff systems. 


7.3.3 Spliting the computational effort 


The points were the approximate solution must be computed are: 


for yn, 1 point: 
for Yn, 2 points : 
for Yn; J points : 


t; ṣ- 
t; ṣ- 


ti + 


The number of those points is 7(y + 1)/2. 


Proposition 7.1. 


L H 
Z +H 
L o u+ aa‘, GAH 


most İ(g + 1)/2] processors, without communications. 
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The computing effort can be uniformely distributed between at 


Proof. Note that one processor must evaluates the points of the last line. When it is 


not any restriction of the communication, the number of processor must divides the 


value j(j + 1)/2. 


Examples: 


- 


we. 
ll 


we. 
ll 


Processor 0: 
Processor 1: 
Processor 0: 


o> 
. 
| 
T 


~ 
öz 
| 
T 


© 

| 

T 
azem] əlim 


Processor 1: £, + 
Processor 1: £, + 
Processor 0: t; + 
Proc. 0: t+ H 
Proc. 1: t+ a ; 
Proc. 1: #+4 , 
Proc. 0: t; +Ë , 
Proc. 2: ti + a ; 


tt +H 

? tita 
2 

? t+ > 


aliss s 


7.3.4 Outline the algorithm 


We propose the following algorithm: 


, +H 

, t+ —— 
ti +H 
tit , +H 
ti + H ? t; - AH ? 


Step 0. Find the number of processors, P, according to the requested approximation 
order. Set m = 1. 


Step m. Each processor p (with p = 0,..., P) follows the substeps: 


(a) find the approximations in the points for which is responsable; 


(b) send the computed dates and receive dates from all the other processors; 


(c) find the final solution of the step, according the relation (7.1). 
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(d) m — m+1. Ifm < N, continue with the next step, else stop the algorithm. 
Remarks. 


1. The maximum order obtained with P processors is 2P + 1. The degree of the 
parallelism of the proposed algorithm depends on the requested order of accuracy. 


2. The approximation error produced by the parallel algorithm is the same like of 
the sequential algoritm. Many parallel methods proposed for the integration of 
differential systems does not satisfy this request. 


3. The value T;; can be express as a linear combination of the yn, values: 


.— | L.. 
Tij = Yn + C2Yn2 J F Cpln)- 


From the order conditions it results a linear system from which we get the weight 


values: 
cı Fey ede, el, 
C2 4 4 Cj 
€ + — ... — = 
1 4 İB ? 
| C2 | Ci — 
Cı Ak-2 İ İ (7) üb 
cı beş 4 + (—1)’e; = 0 


The last equation is the result of the stability request at infinity. For example, 


13 1 45 


~ eq Tn + m. 


T33 — ? 


Expressing T}; as a linear combination of some known values, the time necessary 
for the stage (c) is reduced comparing with the recursive evaluation time. 


7.3.5 Theoretical study of the efficiency 
. Analysing the algorithm, we note that: 


e The efficency of the parallel implementation of the algorithm depends on the 
ratio (Ta + 7,)/T-, where T, is the time necessary to acomplish the phase x 
(x = a,b,c) at some step. 


e The numerical efficiency depends on the ratio between the computing times for 
the ya, (to + H) values and the time for computing the final extrapolated value 


Un. 
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7.4 Second algorithm 


7.4.1 Problem to be solved 


We consider a stiff system of the form 
y(t) = A()y(t) + g(t), 1 e€ [0,7], 


with the initial condition y(0) = y°. Note with D the dimension of this system. The 
class of such stiff systems is not empty (see (S17) and (518) from Appendix 1). For 
the numerical integration, the interval is divided in N equals subintervals of h lenght. 


7.4.2 "The basic method 
Applying the implicit Euler rule 
Ynti = Yn thfns, n=0,...,N—1, 


we get 
[I m hA(tn41)|Yn4a = Yn + hg(ta+1). 


One step n of the sequential algorithm consists on the following: 
1. Evaluate A(tn41) and g(tn41). 


2. Solve the linear system. If we use a Gauss like procedure, the principal phases 
in solving a linear system Qg = 6 are: 


(a) Transform Q to an upper triangular form. 
(b) Transform b in the same manner. 


(c) Solve the upper triangular linear system in z. 


Remarks. The phases (1) and (2a) are independently on the initial value. 


7.4.3 Splitting the computational effort 


We dispose of a system with P processors, connected in a circular network. We choose 
two integer values K and r so that 


KPr — N — Tİh, 


and we make the distribution of the points of the integration interval as in Figure 7.1. 
The idea is to compute at each £ stage, (1) and (2a) in parallel, and (2b) and (2c) 
in serial mode. 
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Figure 7.1: Responsability of the processors 


7.4.4 Outline of the algorithm 


There are K stages. In stage k, where k = 0,...,K — 1, the processor p (p = 
0,..., P — 1) must execute the following: 


1. for ? = 1,...,r 


(a) evaluate A(tkrPh+prh+jh) and g(tkrPh+prh+jh); 


(b) find the matrix C; so that C;[I — RA(tkrPh+prh4;n)| has a triangular upper 
form; 


2. if p #0 or p = 0 and k F 0, receive from processor (p — 1) mod P the vector 
YkrP+pr, else the entry vector is yo; 


3. for ? = 1,...,r 


(a) transform the input vector applying C;; 


(b) find the solution of the transformed system, YprP+pr+i3 


4. ifp#P-lorp=hK-landkF K —1, send ypypy(p41)r to the processor p+ 1 
and set k — k + 1, else stop the algorithm. 


7.4.5 Theoretical study of the efficiency 
Note: 
F - the medium time necessary to evaluate A(-) and g(-); 
U - the medium time necessary to determine a transforming matrix C; 


V - the medium time necessary to transform a D-dimensional vector according to the 
matrix C and to find the solution of the transformed upper triangular system; 
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$ - the medium time necessary to send a D-dimensional vector to a processor directly 
linked with the current processor; 


R - the medium time necessary to receive a D-dimensional vector from a processor 
directly linked with the current processor; 


G - the medium time necessary to solve a linear system of dimension D with the 
standard GAUSS procedure. 


The time evolution of the first algorithm stage looks like in Figure 7.2. With a continue 
line we have plotted the time when a processor is busy. 

The condition that the processor 0 does not wait until the processor P — 1 finishis 
the first stage is 


rU+F)4rV4+S4+rU+Fh)>r(U+P)+(P-D0rV45)4+(P-2)R+R4rV4+8 
that means ULF 
TV” 
(S + R)(P -1) 
~U+F-(P-1)V’ 


If these inequalities are satisfied, then there are not waiting times in the stage k > 0 


(Figure 7.3). 


or 


r 


U+F>(P-DV. 


Stage 0 
0 r(U +F) V S Stage 1 


r(U SE F) ? IR dbi Stage 1 


P— r(U + F) k Stage 1 
Pai dt 


Processor 


Figure 7.2: The first stage, in time 


Theoretically, we have the following equations [111]: 


T,=A[r(U+ FP)4rV4+S54+ R]4+(P-D0V 4 S) 4 (P — 2)9, 
Tı = PKİV(U + F)+rV]= PKr(U-—F-V), 


and then 

Ti 1 

E K(R+S)+(P-1)(rV 4+ S) EP — 2)57 
KrU+F+V) 
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Figure 7.3: The stage k > 0, in time 


Note that 
Tı =To+N(U+V—-G). 
Remarks. 
1, S x R and S does not depend on D; 


2. U, F and V depend superlinear on the system dimension. In addition, F depends 
on the system function complexity. 


3. The stages number K has influence on the efficiency. For system with a small 


D, K must be small, because U + F x S. 


4. The approximation error produced by the parallel algorithm is equal to the one 
produced by the sequential algorithm. 


5. To maximize E5"”, for a given system (for which we know the values 5, T, U, F), 
and a given number of integration steps, N, we search (P, K,r) for which 


K(S + R)+(P-1)\(rV ES) (P-2)R 
(K,Pər)€A KrU+F+4+V) 


pi 


where 
A — ((K, P,r) € N? [rU +F —(P — Vİ > (S+ R)(P-1), PKr= N} 


For a gived P value, the optimal parameters are the followings: 


S+R — 1 a 
Tmin — m O N ——, K min m əə. N — 
V VP(P-—1) S+R P 


6. For simple system function, F € U, F € V, and U/V < (D + 1)/4. The 
condition that the processor 0 has not waiting times implies that 
D+3 


< —. 
pro 


Chapter 8 
NUMERICAL RESULTS 


The purpose of this chapter is to emphasize the implementation performances of the 
methods proposed in the last chapters. 

The efficiency of a method can be measured in many ways. Between the computer 
independent measurements, we mention the number of function evaluations, the num- 
ber of JACOBIAN’s matrix evaluations and the number of matrix inversions. Between 
the computer dependent mesurements we mention the execution time of some codes 
associated with the methods. 


8.1 Split ADAMS-MOULTON schemes 


In this section we analyse the followings: 


e the performances of the split ADAMAS-MOULTON scheme of second order, com- 
pared with those of the classical ADAMS-MOULTON’s formula with two steps, 
which supposes a reduced computing effort, but it is unstable for reasonable 
stepsizes in the case of the integration of a stiff system; 


e the performances of the split scheme with three steps, compared with the split 
backward differentiation formula with three steps, which supposse the same com- 
putational effort. 


Similar tests are presented in [99], [112] and [114] (the test systems for which the 
integration with the split schemes is succesfully are (96), (S17) and (521) from the 
Appendix 1). The split methods are compared with an extrapolation method. 


8.1.1 Experiment 1. Comparison between the split scheme 
and the classical formula 


The test systems are the followings: 


1, the linear system (52) from the Appendix 1 integrated on the interval [0,100], 
with the stepsize h = 0.5, and starting from y(0) = (1,0)7. The ratio of the 
extreme eigenvalues is S = 10°; 
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2. the nonlinear system ($20) from the Appendix 1 integrated on the interval [0,20], 
with the stepsize h = 0.2, and starting from y(0) = (1,1,1,1)*. 
The numerical methods are the followings: 


AM?: the ADAMS-MOULTON’s formula with two steps. The starting value for the 
NEWTON’s iterations for solving the implicit equation is given by the ADAMS- 
BASHFORTH s explicit two step formula; 


SAM: the split ADAMS-MOULTON scheme with two steps (5.7) and with 9 = —1/3 
(for this value the predictor formula is A-stable). The scheme is A-stable. The 
starting value for the NEWTON’s iterations is the same like for the classical 
method. 


EX: for a high approximation of the nonlinear system solution we use the 
DEUHLHARD’s four order extrapolation method based on the linear-implicit EU- 
LER rule (3.21), which is A(89.77°)-stable. 


The starting values for the two steps formulae are the followings: 
1. in the case of the linear system, the exact values of the solution, 


2. in the case of the nonlinear system, the approximate values produced by the 
extrapolation method. 


The NEWTON’s iterations are simplified, i.e. we use only one JACOBIAN’s matrix 
inversion, the evaluated one for the starting value. 
The results are the followings: 


1. for the first system, we draw the conclusions that 


(a) the AM; formula is unstable in both solution components; 


(b) the SAM: scheme successfully integrates the system on the given interval. 


Note the absolute error e = (€1, €2,..., €m)" where 


e(t) = |yj(t) — (yen) sl, te {tilk =1,..., N}, J=l1l,..., M. 
For example, in Figure 8.1 we have plotted the linear interpolation curves of the 
discrete values to sustain the above statements. 
The computational effort was measured. At one step we have count: 
AM?: one JACOBIAN’s matrix evaluation, one matrix inversion and two function 
evaluations; 
SAM): one JACOBIAN’s matrix evaluation, one matrix inversion and three func- 


tion evaluations. 


2. for the second system, we draw the conclusions that 
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Figure 8.1: Error behaviour in the second solution component 


(a) the AM? formula is unstable in the last two solution components (stiff 
components); 


(b) the SAM, scheme successfully integrates the system. 


For example, in Table 8.1 we present some discrete values of the solutions. 


Applying the split scheme, for a precision € = 1071? of the NEWTON’s simplified 
iterations, there are necessary at each step: one JACOBIAN’s matrix evaluation, 
one matrix inversion and nine function evaluations (in mean). 


8.1.2 Experiment 2. Comparison between the split ADAMS- 
MOULTON scheme and a split backward differentiation 
scheme 


The test systems are the followings: 


1, the linear system (81) from the Appendix 1 integrated on the interval [0,20], 
with the stepsize h = 0.05, and starting from y(0) = (1,1,1,1)7. The ratio of 
the extreme eigenvalues is S = 2 - 107. 


2. the nonlinear system (522) from the Appendix | integrated on the interval [0,0.6], 
with the stepsize h = 0.002, and starting from y(0) = (1, 0, 0)7. The JACOBIAN’s 
matrix eigenvalues depend on t: A,(0) = A2(0) = 0, As(t) = —0.04, A1(40) = 
0, A2(40) = —3100, \3(40) = —0.03: 


3. the nonlinear system ($29) from the Appendix 1 integrated on the interval 
[0,150], with the stepsize h = 0.5, and starting from y(0) = (1,1,0)7. The 
stiff ratio depends on t. For the starting values S(0) = 1.375 - 10”. 


The numerical methods are the followings: 
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Table 8.1: Approximate solutions of the nonlinear system 


TUR 


1.856903 
0.033898 
0.034356 
0.034451 
1.999940 
0.0399973 
0.040014 
0.040030 
2 
0.04 
0.040016 
0.040032 


1.864032 
0.034188 
0.0344793 
0.033156 
1.999955 
0.039998 
0.040014 
0.040030 
2 
0.04 
0.040016 
0.040032 


1.863837 
0.034183 
0.373462 
3.131642 
1.999954 
0.039998 
54.65505 
7682697.13 
2 
0.04 
31365.890 
7.568 - 10'4 


2110 
101 50 
20 1100 


SAMa: the split ADAMS-MOULTON scheme with three steps and of four order (5.9) 
with 9 = —7/24, in which case the scheme is A(89.296°)-stable. The star- 
ting value for the NEWTON’s simplified iterations are produced by the ADAMS- 
BASHFORTH s three-step explicit formula, 


SBDFs: the CASH”s split backward differentiation scheme (2.35) with three steps and 
third order (which supposes a similar effort like the above method) for 6 = —1., 
in which case the scheme is L-stable. The starting value for the NEWTON’s 
simplified iterations are produced by the EULER’s explicit rule; 


RK: the standard RUNGE-KUTTA’s process of four order, implemented for finding 
the stiff components of the systems; 


EX: the analougous extrapolation method to that of the first experiment, but of five 
order, for a high approximation of the exact solution for comparing the above 
methods. 


The starting values for the two steps formulae are the followings: 
1. in the case of the linear system, the exact values of the solution; 


2. in the case of the nonlinear systems, the approximate values produced by the 
extrapolation method. 


Like in the first experiment, the NEWTON’s procedure is implemented in a simpli- 
fied version with only one matrix inversion. The tolerated error between two succesive 
iterations is e = 1077. 

The results are the followings: 
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Figure 8.2: Absolute errors in the first system stiff components 


1. for the first system, we draw the conclusions that 
(a) the RK process can not integrate the stiff components, y3, y4, though it 
gives satisfactory results in the first two components; 
b) the split schemes successfully integrate the system; 
P y 8 y ; 


c) the error produce e approximation of the exact solution wi e 
th produced by the approximati f th t soluti ith th 
SAMs scheme is lowest than the one produced by the SBDF3; 


(d) the computational effort is the same for both split schemes: three function 
evaluations, one JACOBIAN’s matrix evaluation, and one matrix inversion. 


For example, in Figure 8.2 we have plotted the absolute error curves. 


2. integrating the second system, the RK process gives an unstable solution (ap- 
proximate values for only 12 steps). The stiff component is yə. In Figure 8.3 we 
have plotted the approximate solutions produced applying the two split schemes 
and the absolute error curves. 


3. for the third system, we draw the conclusions that: 
(a) the solution produced by the RK process is unstable in all solution compo- 
nents; 
(b) the error produced by the split CASH”s scheme is unbounded; 


(c) the split ADAMS-MOULTON scheme successfully integrates the system and 
the mean computing effort at one step can be express in nine function 
evaluations, one matrix inversion, and one JACOBIAN’s matrix evaluation. 


For example, in Table 8.2, we present some significant values. 
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Figure 8.3: Approximate solutions and absolute error in the second component of the 
exact solution 


For the instability of the other schemes, we can mention some values: 


7668.407 279679399.97 
ypo — | 329.1616 | , yBPFss — | 11947028.66 | , 
423.4813 15687767.34 
—3673 7.044 - 109€ 
yil = | 6.413 İ, iln = | 1.525 - 1040 
6.936 1.943 - 1040 


8.1.3 Conclusions 


The split ADAMS-MOULTON formulae can be used with success in the integration of 
stiff systems. Altrough using a such scheme is more expensive in the sense of the 
computing effort, the advantage is the possibility to integrate some systems for which 
the classical ADAMS-MOULTON’s formulae reach some stability difficulties. 

Comparing the split ADAMS-MOULTON scheme with the CASH’s split backward 
differentiation scheme with the same stepnumber, one can see that the first one pro- 
duces a discrete solution more closed to the exact solution than that produced by the 
second one. 

The split scheme can be implemented in a numerical code which: 


1. uses an ADAMS-MOULTON’s formula as basic formula; 


2. includes an instability test which determines the situation when we deal with a 
stiff system; 


3. changes the basic formula with the split corresponding scheme when the insta- 
bility test is successfully. 
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Table 8.2: Approximate solution produced by SAM integrating the third system 


vüsal Tu" [AOS ayi 


1.251492 1.254754 0.003262 
5 | 10 1.070882 1.070882 6.164 - 1077 
0.608467 0.608069 0.000399 
1.625682 1.626610 0.000927 
50 1100 1.558227 1.559085 0.000858 
7.302124 7.300877 0.001247 
0.819860 0.820188 0.000328 
150 | 300 0.950671 0.951237 0.000565 
20.38286 20.39756 0.014700 


A such implementation version has an important advantage: the stored values from 


one step to another must be not modified when we change the formula. 


8.2 Split second derivative scheme 


In this section we analyse the followings: 


e the performances of the split ENRIGHT scheme with three steps (5.17), compared 
with those of the classical ENRIGHT’s second derivative formula with three steps 
(5.2), which has some difficulties to integrate the stiff systems with some eigen- 
values closed to the imaginary axis; 


e the performances of the split second derivative backward differentiation scheme 
with four steps (5.20) compared with those of the classical second derivative 
backward differentiation formula with the same stepnumber (3.3) which is only 
stiffly stable. 


For the chosen @ values, the split schemes are A-stable. 
Similar tests are presented in [109]. 


8.2.1 Experiment 1. Classical formula instability 
The test system is a modification of the system (89) from the Appendix 1: 
{ yt) = ayı (t) + byəit), 
yət) = —byr(t) + ayr(t) 


The initial values are y(0) = (0,1), and the integration interval [0,7]. The exact 
solution is y(t) = (e“f sin bt, e” cos bt). 


We propose the following values: 
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1. for the split ENRIGHT scheme, h = 0.5, a = —0.04, b = 3.6, T = 50; 


2. for the split second derivative backward differentiation scheme, h = 0.5, a = 


—0.05, b= 12, T = 10. 
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Figure 8.4: Exact and approximate solutions produced by EN; and SES; in the first 
component 


The numerical methods are the folowings: 


ENs: the ENRIGHT’s second derivative formula with three steps and of five order 
(5.2). The starting value for the simplified NEWTON’s iterations, for solving the 
implicit predictor formula, is the one produced by the EULER’s explicit rule. 


SES3: the split ENRIGHT scheme with three steps and of five order (5.17) with 
0 = 1/12 (for this value the scheme is A-stable). The starting value is the same 
like for the above mentionated method. 


SBDEFA: the second derivative backward differentiation formula with four steps and 


of five order (5.3). 


55BD54: the split second derivative backward differentiation scheme, with four steps 
and of five order, (5.20), with @ = —1/3 (for this value the scheme is A-stable). 


As starting values for the multistep formulae we have take the exact values of the 
solution. 

The NEWTON’s procedure for solving the implicit equations is built-up in a simpli- 
fied variant with only one matrix inversion (the JACOBIAN’s matrix evaluated at the 
starting value). 

We have get the following results: 


1. the EN; and SBDF, formulae are unstable in the both solution components; 


2. the SES; and 55BD5, schemes successfully integrate the system. 
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Figure 8.5: Exact and approximate solutions produced by SBDEF, and SSBDSy, in the 
second component 


For example, in Figures 8.4 and 8.5 we have plotted the exact solution and the ap- 


proximate ones. 
The computational effort to one step is the following: 


- for EN; and SBDF4, one JACOBIAN’s matrix evaluation, one matrix inversion, two 
function evaluations and one derivative evaluation; 


- for SES; and SSBDS4, one JACOBIAN’s matrix evaluation, one matrix inversion, 
three function evaluations and two derivative evaluations. 


8.2.2 Experiment 2. Error of the split schemes 
The test systems are the followings: 


1, the linear system (82) from the Appendix 1 integrated on the interval [0,20], 
with the stepsize h = 0.5, and starting from y(0) = (1, 0)7, 


2. the nonlinear system ($29) from the Appendix 1 integrated on the interval 
[0,100], with the stepsize h = 0.5, and starting from y(0) = (1, 1,0)". 


The numerical methods are the same like in the first experiment. 
We have get the following results: 


1, the errors (absolute or relative) produced by the approximation of the exact 
solution with the split schemes are lowest than the ones produced by the classical 


formulae; 
2. the relative error produced by 55BD5A is lowest than the one of the SESs. 


To prove these statements, in Figure 8.6 we have plotted the absolute error curves 
for the first system, and in Figure 8.7, the relative error curves for the second system. 
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Figure 8.6: Absolute errors approximating the first system solution 
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Figure 8.7: Relative errors approximating the second system solution 


8.2.3 Conclusions 


The results of the first experiment confirm that the stability domains of the split 
schemes are largest than the ones of the classical corresponding formulae. 

The split schemes produce some approximation errors lowest than the ones cor- 
responding to the classical formulae. Comparing the split ENRIGHT scheme with the 
split second derivative backward differentiation scheme with the same order (but with 
different stepnumbers), we can see that, for the second one, the error has some smallest 
discrete values. 

The proposed split schemes are proving their efficiency in some implementations 
like the following: 


1. use a second derivative multistep formula as basic formula; 


2. include an instability test which identifies the situation when the stiff system 
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does not correspond to the stability domain of the basic formula; 


3. switch from the basic formula to the split corresponding scheme if the instability 
test is successfully. 


8.3 Exponentialy fitted formulae 


In this section we study the performances of the parameter dependent backward diffe- 
rentiation formula with two steps, compared with two GEAR’s backward differentiation 
formulae, one of the same order, the other of the same stepnumber. 

Applying the method of exponential fitting, we expect that the parameter depen- 
dent formula produces the smallest error. 

Similar test are presented in [100]. We must mention the very good results in the 
integration of the linear system (512) and the nonlinear system (S23). 


8.3.1 Experiment 1. Backward differentiation formulae ap- 
plied to a nonhomogeneous equation 


The test equation is (S10) from the Appendice 1, which is integrated on the interval 
[0,100], with the stepsize h = 0.1, and starting from y(0) = 10. The exact solution is 
y(t) = E(t) + 106770". 


Note the following numerical methods: 


BDF: the GEAR’s onestep backward differentiation formula (1.12) of first order (the 
EULER’s implicit rule) which is A-stable; 


BDF»: the GEAR’s two-step backward differentiation formula (2.31) of second order 
which is also A-stable. The starting value for the multistep formula is the exact 
value of the solution. 


BDP?: the parameter dependent backward differentiation formula (6.3) with two 
steps and of first order which is also A-stable for an a value chosen by the 
request of exponential fitting: 


24 —3 €? 4+ del —] 
a = a(q) = (74-35) 46 


(ey 
where 
1744771 altfel, 


The switch between the two fitting points is made when the stiff component, 
10e~7°", no more affects the solution: ¢ is chosed so that 10e < 1077”. 


EU: the EULER’s explicit rule for proving the stiff character of the system. 
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We draw the conclusion that, for the test equation, the exponential fitted formula 
produces an approximation error lowest than that produced by the two GEAR’s for- 
mulae. The unstable behaviour of the EU’s solution le to the conclusion that we deal 
with an stiff equation. 


In Figure 8.8 one can see the plots of the absolute error curves. 


0 05 1 15 2 25 3 35 4 005 115 225 3 35 4 
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Figure 8.8: Approximation errors produced by the backward differentiation formulae 


8.3.2 Experiment 2. Backward differentiation formulae ap- 
plied to a perturbated linear system 

We have tested the methods from the first experiment in the integration of the system 

(S13) from the Appendice 1, on the interval [0,25], with the stepsize h = 0.02. The 


eigenvalues ratio is $ = 1.5 - 10°. 
The a parameter is chosed as 


_ (2q — 3)e"4 + 4e1 — 1 
a= a(q) = 6-0 


where q = —h. 

We draw the conclusion that, for the given system, the parameter dependent for- 
mula produces a lowest error than the one of the GEAR’s formulae. The EU formula 
has an unstable behaviour for the given stepsize. For proving this statement we have 
plotted in Figure 8.9 the absolute error curves in each solution component. 


8.3.3 Experiment 3. Backward differentiation formulae ap- 
plied to a nonlinear system 
We draw the same conclusion like for the last two experiments, when we integration the 


nonlinear system (520) from the Appendice 1, in the interval [0,20], with the stepsize 
h = 0.2, and the starting values y(0) = (1,1,1,1)7. 
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Figure 8.9: Absolute errors in the approximation of the perturbated linear system 
solution 


The method which give us an high approximation solution is the DEUFLHARD’s 
extrapolation method of four order, based on the linear-implicit EULER rule (3.21) 
(noted EX). We use also the standard RUNGE-KUTTA’s process of four order for 
finding the stiff components of the system (noted RK). 

The a parameter is chosed as 


29 4 q 
a x atq) — 174 3)67 + de L 
(eT) 
where q = —100 - 2. 

The RK method application indicates that the third and the fourth components of 
the exact solution are stiff components. 

For proving the above conclusion we have plotted in the Figure 8.10 the absolute 
errors in two solution components with disctinct character (the first one is a smooth 
component, the second one, a stiff component). 

The computing effort is the same for both two step formulae. 


8.3.4 Conclusions 


From the point of view of the aproximation error, an exponential fitted backward 
differentiation formula is more proper for the integration of a stiff system than the 
classical GEAR”s formulae. 


8.4 Nonlinear formulae 


We study the efficiency in the implementation of one exponentially fitted nonlinear 
method presented in Chapter 5, compared with closed formulae: 
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Figure 8.10: Absolute errors in the approximation of two distinct components 


ə the linear-implicit EULER’s rule; 
e the linear-implicit midpoint rule; 
e the EULER’s explicit rule. 


The comparison between some A-stable nonlinear explicit formulae and some 
stiffly-stable linear multistep formulae can be done only from the error point of view. 
Comparing the computational effort, we draw the conclusion that the explicit nonlinear 
formulae are more convenient. 

Similar tests are presented in [103] and [115]. The discussed methods are those with 
exponentially fitted coeficients, and the comparing methods are the trapezoidal rule, 
the GEAR’s backward differentiation formulae and the ENRIGHT’s second derivative 
formulae. We can mention the improvement in error decrease using the nonlinear 
methods in the integration of the system ($14) and (S22). 


8.4.1 Experiment 1. Integrating a separable system 


We consider the test system (S21) from the Appendix 1 on the integration interval 
10,201, with the stepsize h = 0.1. The initial ratio of the JACOBIAN’s matrix eigenvalues 
is S = 106. 


Note the following numerical methods: 


NL: the nonlinear method with two steps (6.18), with maximum third order, which 
is A-stable. The starting value is the exact solution value at £ = h; 


EL: the linear-implicit EULER’s rule (3.21) of first order, which is A-stable; 


EE: the EULER’s explicit rule of first order, for finding the stiff components of the 
exact solution. 
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We have get the following results: 


1. the behaviour of the EE solution lies us to the conclusion that the first two 
components of the exact solution are stiff components; 


2. the nonlinear formula produces the exact solution. 


For proving these statements, we present in Table 8.3 some values of the numerical 


solutions. 


Table 8.3: Exact and approximate solutions of the separable system 


—1.077 - 107 — i 
—1.675 - 1074 
—7.75 
—0.8455 
—1.025 - 10714 


—4.803 - 10714 
—9.93 
—0. nan 


0 
—10 
—0.04783 


From the calculus effort point of view, we mention that, at each step, there are 


necessary the followings: 


1. for NL, one JACOBIAN’s matrix and one function evaluations; 
2. for EL, one JACOBIAN’s matrix and one function evaluations; 


3. for EE, one function evaluation. 


8.4.2 Experiment 2. Integrating a nonlinear system 


The test system is (820) from the Appendix 1, integrated on the interval [0,20] with the 
starting values y(0) = (1, 1, 1, 1)7, and the stepsize h = 0.2. Note that the JACOBIAN’s 
matrix eigenvalues are the same in all the integration interval. The ratio between 
the largest and the smallest one is S = 107, The JACOBIAN’s matrix is diagonally 
dominated. 


The numerical methods are the followings: 


NL: the above mentionated nonlinear formula. The starting value is produced by the 


EX method; 
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EX: the DEUFLHARD’s extrapolation method of four order for a high approximation 
of the exact solution; 


PM: the linear-implicit midpoint rule, which is A-stable, and has second order; 

EE: the EULER’s explicit rule, for finding the stiff components of the exact solution. 
We draw the following conclusions: 
1. the EE solution is unstable in the last three components (stiff components); 
2. the nonlinear formula produces the exact solution in the first component; 


3. the approximation error produced applying the NL formula is lowest than the 
one produced by PM, in the first half of the integration interval. 


The computational effort is evaluated at one step: 


1. for PM are necessary one function evaluation and two multiplication between a 
matrix and a vector; 


2. for NL are necessary one JACOBIAN’s matrix evaluation and one multiplication 
between a matrix and a vector. 


To prove the above statements, we mention in Table 8.4 the approximate values of 
the exact solution. 


Table 8.4: Exact and approximate solutions of the nonlinear system 


2 


0.0400 . 
u(1002) = yEX = 0.4002 l Yio = 
0.0416 
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8.4.3 Conclusions 


The tested nonlinear formula can successfully integrates some stiff systems. The ap- 
proximation error is comparable with those produced by some implicit A-stable linear 
multistep formulae. 


8.5 Onestep hybrid methods 


In this section we study the performances of two hybrid methods proposed in the 
Chapter 6, compared with an ENGLAND’s method [37] with the same computational 
effort. 

We mention that in the papers [98], [102] and [106] are given the test resulting from 
other three hybrid methods proposed in Chapter 6. The compared methods are the 
ENRIGHT’s onestep second derivative formula, the OBRECHKOFF’s second derivative 
formula, and a BOCKHOVEN’s hybrid scheme. The test systems are (512), (513), (S17) 
and ($25) from the Appendice 1. 


8.5.1 Experiment 1. Integrating one stiff equation 


The test equation is (S10) from the Appendice 1, integrated on the interval [0,20], 
with the stepsize h = 0.5, and the starting value y(0) = 10. The exact solution is 
y(t) = E(t) + 106770". 


2 4 6 8 10 12 14 16 18 20 13 14 15 16 17 18 19 20 
t t 


Figure 8.11: Relative errors in the integration of the stiff equation 
The tested numerical methods are the followings: 
H1: the ENGLAND’s L-stabil scheme of third order (6.19); 
H2: the L-stable hybrid method of third order (6.25); 
H3: the strong A-stable hybrid method of third order (6.24); 
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RK: the standard RUNGE-KUTTA’s process of four order. 


The starting value for the NEWTON’s procedure for solving the implicit equations 
is given by the EULER’s explicit rule. 

We draw the conclusion that the first three methods successfully integrate the 
system. The lowest error is produced by H2 in a first integration interval, and by H3 
in the second one. In Figure 8.11 we have plotted the relative errors curves. 

The computational effort is the same for the first three methods and consists, at 
each step, in two function evaluations, one JACOBIAN’s matrix evaluation, and one 
matrix inversion. 


8.5.2 Experiment 2. Integrating an linear system 


We test the same numerical methods on the system (52) from the Appendice 1, in- 
tegrated on the interval [0,20], with the stepsize h = 0.5, and the starting values 
y(0) = (1,0)7. The exact solution of the problem is xy, (£) = 2e7f — e712, yo(t) = 
—e”t ze” 10008, 


We get the following results: 


1. the H1 and H2 methods produce approximately the same absolute error (the 
small difference between them indicates that the H2 solution is more closed to 
the exact solution); 


2. the H3 method produces a lowest error than the ones of the H1 and H2 methods 
only in a second integration subinterval; 


3. the RK can not integrate the problem. 


In Figure 8.12 one can see the plots of the relative errors curves. 
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Figure 8.12: Relative errors of the approximative solutions 
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Table 8.5: Approximative solution in the integration of the nonlinear system 


in — mill vite 
1.181287 0.004596 1.181286 0.004596 1.181 
0.122202 0.096288 0.122204 0.096286 0.3427 
0.058113 0.090867 0.061166 0.087815 85.62 
—0.055171 | 0.064238 1 —0.054026 | 0.063094 5515 
1.632159 0.012007 1.632159 0.012007 1.632 
0.025358 0.000739 0.025359 0.000739 0.02946 
0.263383 0.011295 0.263386 0.011298 | 1.267 - 1010 
0.027210 0.001211 0.027210 0.001211 | 1.078 - 107? 
1.864693 0.008383 1.864693 0.008383 1.865 
0.034215 0.000694 0.034215 0.000694 0.03424 
0.346532 0.006399 0.346532 0.006399 | 2.071 - 107° 
0.035929 0.000699 0.035929 0.000699 | 1.611 - 1099 


Since the system is linear, the computing effort is the same for the first three 
methods. At one step, there are necessary two function evaluations, one JACOBIAN’s 
matrix evaluation, and one matrix inversion. 


8.5.3 Experiment 3. Integrating a nonlinear system 


We test only ET and H2 schemes in the integration of a nonlinear system, in the aim 
to prove that the produced errors are the same or very closed each another. The test 
nonlinear system is (S20) from the Appendice 1, integrated on the interval [0,2], with 
the stepsize h = 0.2. 

For a high approximative solution we use the DEUHLHARD’s extrapolation method 
of four order. 

The starting value for the NEWTON’s procedure for solving the implicit equations 
is that produced by the EULER’s explicit rule. 

The results confirm the above statement (see Table 8.5). 

There are no reported differences between H1 and H2 scheme in the computational 
volume. 


8.5.4 Conclusions 


We have numerical proved that the first third order hybrid method, proposed in Chap- 
ter 6, is comparable with the ENGLAND’s hybrid method of the same order in both 
senses, of approximation error and computational effort. 

The second scheme is proper for integration when we are intersted only in the 
solution at the final integration value. 
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8.6 Parallel algorithms 


In this section we study the performances of the parallel algorithms proposed in Chap- 
ter 7. The tests was made on a transputer system (T-800) under the PARIX operation 
system. The results of these tests are also presented in the paper [110]. 


8.6.1 Experiment 1. The efficiency of the parallel extrapo- 
lation algorithm 


The test systems are the followings: 
l. a generalization of (S6) from the Appendix 1: 
yi(t)=—Py(t), te l011 t@=1,...,D. 
For D = 10 we get (S6). 
2. the system ($40) from the Appendix 1. 
The implemented algorithms are the followings: 
e the parallel algorithm from section 7.3.4; 


e the sequential DEUFLHARD’s extrapolation method (with minimum number of 
arithmetical operations) 


İp 


Figure 8.13: The running time of some algorithm implementations (a) for P = 2 
processors (b) for P = 3 processors 


We study the effect of changing: 


1. the system dimension D; 
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2. the stepnumber N; 


3. the processor number P. 


We have measured the folowings times: 


- the execution time of the sequential code, To; 


- the execution time of the parallel code using P processors, Tp: 


- the medium computing time of a processor, T$. 


We present the results in a graphical form. 


(a) 


0 2 4 6 8 10 12 14 16 18 20 
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Figure 8.14: The efficiency E(D, N) for some (a) step numbers (b) dimensions 


Analysing the figures we draw the following conclusions: 


1. 


By the implementation of the parallel algorithms, the running time of the extra- 
polation method can be improved (Figure 8.13). 


. The efficiency Eş of the parallel algorithm, which use P = 2 processors, depends 


on the dimension D of the system (Figure 8.14.a), but not on the stepnumber 
N (Figure 8.14.b). For D > 10, E; is closed to the ideal value; 


. For some given values N and D, the efficiency has not a significant dependence 


on the number of processors, P (Figure 8.15.a). This fact is not valid for most 
parallel algorithms. The explanation lies, in our case, to the fact that the com- 
puting time of the sequential algorithm, 7o, depends on the requested accuracy 
of the numerical solution, which is a function on P. 


. For some given value N and P, the total computing time increases with the 


dimension D so that it covers the negative effects of the communications (Figure 


8.15.b). 
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5. The employing factor U does not depends on the number of processors P and 
on the stepnumber N (Figure 8.16). 


The efficiency, for fixed stepnumber, dimension and processors, depends on the 
complexity of the system function. For instance, for the second system, the efficiency 
of the parallel algorithm, with N = 20000 and P = 2 processors, is 85.191%. 


6—O---- 
— — 49-10V147 _ 
~~ 720 
-1 -0.8 -0.6 -0.4 -0.2 0 
Re(z) 


Figure 8.15: (a) The efficeney Ep(D) for N = 1000 subintervals (b) The communica- 
tions influence for N = 1000 and P = 3 
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Figure 8.16: The employing factor U(D,P) (a) for N = 1000 subintervals (b) for 


P = 2 processors 
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8.6.2 Experiment 2. Parallel algorithm for linear stiff sys- 
tems 


We have tested the proposed algorithm on the system (S17) and on the following 
system 
y(t) = —2°y;(t), te[0,lj, t=1,...,D. 


This system generalizes in D dimension the stiff system (56). 
The implemented algorithms are the followings: 


e the parallel algorithm from section 7.4.4; 


e the EULER’s implicit rule with the GAUSS’s procedure for solving the linear 
resulting system; 


Like in the first experiment, we test the influence on the algorithm performances 
of the N, D and P parameters and of the system function complexity. 

Analysing Table 8.6, we draw the conclusion that the efficiency of the parallel 
implementation of the algorithm increases with the dimension D and decreases with 
the processors number P. For P = 2, EZ” is very closed to the ideal value. 


Table 8.6: Efficiency in the integration of first system with N = 1000 


From Table 8.7 we can deduce that all the measurements of the efficiency increases 
with the number N of the points where the solution of the ODE system is numerical 
evaluated. 

The complexity of the system function has a significant influence on the numerical 
efficiency (compare the values for P = 2 from both tables). For the test systems, the 
value of Enum is closed to 0.8. 

The numerical efficiency is the principal value which has a great influence on the 
efficiency of the parallel algorithm. Hence, Ep depends on D, P and on the complexity 
of the system function. 
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Table 8.7: The efficiency in the integration of the Iserles’s system (2) with the dimen- 
sion D = 2 and P = 2 processors 


8.6.3 Conclusions 


The first proposed parallel algorithm for distributed memory multiprocessors is de- 
signed to solve stiff systems of many ordinary differential equations. It has the advan- 
tage that the approximation error is the same like the one produced by the sequential 
algorithm. 

The second proposed parallel algorithm is designed to solve linear stiff systems with 
variable coefficients and a perturbation depending on time. The approximation error 
is the same like in the sequential algorithm. The efficency of the parallel algorithm 
increases with the dimension and the complexity of the system function and decreases 
with the number of used processors. 


Appendix A 


Classification of the stiff systems 


The numerical methods designed for solving stiff problems are not universally. Some 
numerical codes work better for a particular stiff system subclass, others for another 
subclasses. 

The most popular classification of the stiff systems is due to ENRIGHT et al. [40]: 


Class A: the linear systems with real matrix eigenvalues; 
Class B: the linear systems with complex matrix eigenvalues; 


Class C: the systems with nonlinear coupling between the smooth and stiff solution 
components; 


Class D: the nonlinear systems with a real JACOBIAN’s matrix eigenvalues; 
Class E: the nonlinear systems with a complex JACOBIAN’s matrix eigenevalues. 


For proving the efficiency of a new method designed for stiff problems, the method 
must integrate with success all the stiff systems from STIFF-DETEST package [40]. 

In the next Tables we present the classical stiff systems. Note S the ratio between 
the maximum and minimum absolute eigenvalues evaluated at the starting point. The 
second notation correspond to that of STIFF-DETEST package. 
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Table A.1. Stiff linear systems from Class A 


y! = —10°y, + 0.075 ys 

üş = 7500y; — 0.075y2 

yl = — 1800) + 90042 

yh = izi — 2yi + Yin, 2 2(1)8 

ya = 1000ys — 2000y9 + 10000 

t 1002 — 10y3 + y4 
10y3 — 10y4 


+ koy 

+ kaya 

ya = kiyi + kayo — (kə + ks) ys 

ky = 8.43. 10710, kə = 2.9. 1011, 
ks — 2.46 ı 1019 k4 = 8.76 ı 1076 


zi 
e 
E 


= —100y1 — y2 
—100y3 + ya 
—10000y3 — 100y4 
—10% + aye 

= —ay; — 10y2 

= —4ys 

= —lA 
—0.1y6 


ho 
© 


a 


a = 3.8, 25, 100 | 20 


Table A.3. Stiff perturbated linear systems from Classes A and B 


TO) 


— — I 


= aN — (104 t)e* 


rr er TET 
—6y1 + Əyə + 2sint 
S12 100 10° 
mir. f= yy mp pe 
y! = —4498y; — 5996y2 + 0.006 — t 
y} = 2248.5y, + 2997y2 — 0.503 + 3t 


— UBU; FEUD. Uii — — 
ig = ži FJ 
—100 1000 0 
B= 1 —1000 —100 0 „D = 
0 0 —0.1 


3? 


y’ = UBU"y + Uf ), 
uu = —4, u Uij = Ei 7 J 
1 0 0 
-1 0 0 0 
0 0 —100 —900 1” 
0 0 900 —100 
f(t) = (4? 4:21, 0? — 21, 
800£ + 1, —1000¢ 


5 hu 1 
-20 + 40+ 077 
yi = 0.2(y2 — vi) . 102 
y! = 10y, — (60 — 0.125t)ya + 0.1254 (709 3-53 10 
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Table A.5. Stiff nonlinear systems from Class C 


40041 — 100y2y3 — 300043 
= 302 

—0.04y; + 10*yoy3 
= 0.04y, — 10*y2y3 — 3 - 107v2 


= y3 — 100y1y2 

= ys + 2y4 — 100yi yə — 2 - 108 ys 

= —la" 
va 4 

= —0.013y, — 1000y, y3 
—2500y2 ys 


— 0.01 — (0.01 

(y? + 1001yi + 

— 0.01 — (0.01 

= —yı + 10° y3(1 — y1) 

= —10y2 + 3-107 y3(1 — y2) 
=y = 2 


Table A.7. Stiff nonlinear systems from Class E 


yi = yə 

yy = 50 — yi)ye — Yı 

yi = —(55 + ys)yı + 65y2 
ys = 0.0785(y1 — y2) 


= —7.89 . 107192, — 1.1 . 107y1/3 

= 7.89- 107 yy — 1.13. 10° yoy3 

= 7.89 - 107 yy — 1.1. 107y/ y3+ 1000 
+1.13 - 103), — 1.13 - 105y2y3 


, = 1.1- 107yiys — 1.13 - 103 y4 


1 = l0fyiys + 10*y244 
— 1 — us 
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Table A.8. Some stiff nonlinear systems which model real systems 


No Syste — — — — İ7 1:0” 


10fy + y3 28 yi — Ys 
-ły + yı > Us 

—Ü.01x2 

—y3 yy — ys 

—Y1 — Y3ya 


yə = ((1 — yi) Vİ — vilye — m)/e 
£ — 107 6 


Ü.l,c —1,2 ə, ?—4 
— 250|(—0.0048(ys — 660.2)— 
—0.032( ys — 273.9) — Ly + vəl 
— 0.1(vi — y2) 
2 = 93y1 — 0.26(y3 — va) 
— 0.S7(gə ya) 11(y4 us) 


= —100(y3 — yı) 
— 1/94 [0.9 — 1000(y3 — Ys) — 100y3 (ys — y1)| ; 
= 100(y3 — y1) 
= —100(y5 — us) 
a — 0.99026, 0.99, 0.9 
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Table A.9. Other stiff nonlinear systems which model real systems 


Mofyn 0 — Ti) 


= 77.27İy) + yil — 8.375 - 10 5y; — yə)l 
[ys — (1 + xi )yəl/ 77.27 300 i 
= 0.161(y1 — ys) R. 
ü 


= —1.7ly, + 0.43y2 + 8.32y3 + 7- 1077 
1.71y) — 8.75y2 
= —10.03y3 + 0.43y4 + 0.035ys5 
= 8.32y2 + 1.7ly3 — 1.12y4 
= —1.745y5 + 0.43y6 + 0.43 y7 
= = — 2806s + 0.69y4 + 1.7lys — 0.43y6 + 0.69y7 


7 Ny te! xı £ = 1078 


= =| 
D 
BL — 
Ço 
— 


yy = kiyi + ksys — kayoys — 2k ys 
Y3 = keyiyə — kaya 
y4 = —kayoya 
kı = 1071, kə = 2.9 - 104 
ks — 5 ı 10°, k4 = 104 


? 


= 1.30 i 
j= = 1.80 - 10°[ya — ya(1 + k(n) 
yl, = 1752 — 269ys + 26Ty1 
y4 = 0.1 + 320y2 — 321y4 
= 0.006- e20-7—15000/y1 
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Appendix B 


NEWTON’s iterations error 


B.1 Motivation 


In the numerical treatement of the stiff equations we deal with some implicit formula 
of the following form: 
y—hBfly) = 6 

A convergent method, for which is not necessary that the system function f to saf- 
isfy the inequality |hOf/Oy| < c, is the NEWTON’s method. The classical convergence 
conditions are establish by the KANTOROVICH’s theorem [94]. 

The problems associated with an implementation of the NEWTON’s methods are 
the folowings: 


1. the starting values; 

2. a minimum effort in evaluating the derivative; 

3. a method for solving the associated linear system; 
4. the estimate of the error. 


In this section we are treating only the problem 4. 
Generally, we search the solution of the equation system 


where F: D C R” — R” is a given operator. The NEWTON’s iterations are well 
defined by the process 


att — gf — Fa”y F(a"), k=0,1,... 


where F” is the GATEAUX s derivative. Let ||- || the Euclidian norm on R”. We use 
the same notation for the associated matrix norm. 
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Theorem B.1. (NEWTON-KANTOROVICH’s theorem [94]) Assume that F: D C 
R” — R” is F-differentiable on a convex set Do C D and that 

[E œ) = FY) Sale yll, Ya, y € Do. 


Suppose that there exists an x” € Do such that ||F"(2°)""|| < 9 and a = Byp < §, 
where u > | F'(£°) F(2®)||. Set 
č  (84)7111:- (1 — 2a)H2) t= (89) "+ (= 20)"”7] 
and assume that S(x°,t*) C Do. Then the NEWTON ’s iterates are well-defined, remain 
in S(x°,t*), and converge to a solution a” of F(x) = 0, which is unique in S(x? t) A 
Do. Moreover, the error estimate 
Ie" — all x (9:27 “ay”, k= 0,1... 

holds. 


Note that the theorem assumes that the derivative is LIPSCHITZ continous. 
Counterexample. Let the CAUCHY’s problem 


yi = —-V us, 
y, = -vV Y}, 
yi(1) = yo(1) = 4, 


with at most one solution y(t) = yo(t) = 4/t?, t > 1, which is numerical integrated 
with the EULER’s implicit rule. Assume that, at one step, yn-1 = (4/h?,8/h?)". Note 
Yn = (h? x1, h? x2). Then the system is 


(üzvi 
vir ay =8, 


with the exact solution zı = 4, zə = 0. Therefore y, = (427, 0). The system function 
on the unknown z does not satisfies the LIPSCHITZ’s condition on the region [0,00) x 
İ0, co), but satisfies the HOLDER’s inequality 


|e) - Fy) < yle -= yl’, Yzy € Do, OSpsl 


with y = 3/2 and p = 1/2 [113]. 
In the case of a HOLDER continous derivative, the NEWTON attraction theorem 


1941 tell us only that the iterations converge with an superlinear spead. 
An other theorem for estimating the error is the following. 


Theorem B.2. (NEWTON-MySOVSKII theorem [94]) Suppose that F : D C R > 
R" is F-differentiable on a convex set Do C D and that for each x € Do, Fa) is 
nonsingular and satisfies the LIPSCHITZ condition and İF(ey Tl) < 6,Yz € Do. If 
zü € Do is such that ||F'(£°) "F (2°)|| < u and a = lüyü € 1, as well as S(x°,ro) C 


Do, where 
CO 
J 
ro = 1 ) qt 
j=0 
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then the NEWTON S iterates remain in S(zÜ, ro) and converge to a solution a” of 
F(x) = 0. Moreover, the following estimate holds 


jo” = a*l] < ellet = aP, k= 1,2... 


where 


bu a <ofu(l—o J, k=1,2,... 
j=0 


B.2 Generalization of the classical theorems 


We propose the following generalization of the NEWTON-MysovskII theorem: 


Theorem B.3. /101] Suppose that F: D C R" — R” is F-differentiable on a 
convex set Do C D and that for each x € Do, F'(x) is nonsingular and satisfies the 
HOLDER condition with p < 1 and İF"(ey Vİ < G,Va € Do. If at € Do is such that 


| F’(a°) tt F(2®)|| < u and a = se Oye” <1, as well as S(x°,ro) C Do, where 


~ (p+1)3 -1 
To Ku ) a P 
o 


then the NEVVTON S iterates remain in S(zÜ, ro) and converge to a solution a” of 
F(x) 0. Moreover, the following estimate holds 


je" = 2 || < elie’ — NH, b 1,2... 


where 
oo —1 
Aon, pH)" (p+) -1 7 (pt1)* 
o — — (a PF) 9 7— P ) , k=1,2,... 
pe də 


For the KANTOROVICH’s theorem, we have find a partial generalization. 


Theorem B.4. /97/ Assume that F: D C R" — R” is F-differentiable on a convex 
set Do C D and the HöLDER 5 inequality holds for p < 1. Suppose that there exists an 
x° € Do such that ||F’(x°)7"|| < 6 and a = Bye? < zar, where u > | F’(2°)-* F(2®)]| 
and S(x°, (G-y)71/?) C Do. Then the NEWTON’s iterates are well-defined, remain in 
S(x°, (By)-/"), and converge to a solution a” of F(x) = 0. Moreover, the error 
estimate 


(p+1)} 


1 P 1)? P 
1. — 
P Y P 


IF 


holds. 
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In the case p = 1 we get the NEWTON-KANTOROVICH’s error estimate. 
Example. We apply the last theorem to our above mention counterexample. The 


2 əy (Nth 
71:10 


error estimate is 


3 
for any value x° chosen so that 


E 4 
Mı 


holds [105]. 


Appendix C 


Computing the stability angle 


We propose an algorithm for finding the minimum value of a real function. 


C.1 Motivation 


If the curve, which describes the stability domain boundary of a stiffly-stable formula, 
can be expressed in the parametric form 


0 € (0,27) —> z(9) € C, 


then the stability angle is the value a for which 


where 


Im z(9) 
g(9) = arctg — 5 R z(0) < 0, 
m m z(9) 
5 barışan , Rez(0)50 


Therefore, the problem of computing the stability angle is a problem of onedimen- 
sional minimization. 
Remark. The boundary curve of a linear multistep formula generated by the couple 


(p,a) is . 
ple’) 
o (ei). 


If we deal with a multiderivative linear formula, the characterisc equation has a poly- 
nomial expression and there is a number of curves, equal to the derivative maximum 


0 € (0,27) —> z(9) = 


degree, which can be expressed as some z(@) curves. 
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C.2 Onedimensional minimization 
Let the problem to minimize the real function g : [ao, bo] + R, i.e. finding 


min x). 
r€[ao,bo] g 


When the g expression is very complicated, it is preferable to use a method which 
does not request the derivative evaluation. Such methods are the gold section method 
or the FIBONACCI search method. 

We propose an improvement of those methods. 


C.3 A-algorithm 


The two above mentionated methods can be expressed in a form named A-algorithm 
[94]. For each A € [4,1], the \-algorithm follows the steps: 


Step 1. By a certain method, we determine an interval lap, bo], which includes the 
solution a”. We give the error tolerance e. Set 


Ao = bo — do, Ay — AAo, Az — Ao — A), Yo — do + Ao, Zo = bo — Ao, k =0 
We compute g(yo) and g(zo). 


Step 2. If A, < €, then we put «* = arg min {g(y«), (zx) } and we stop the algorithm, 
else k k+1. 


Step 3. We put Ayo = Ag — Anyi and, if glyk-1) < g(zi-i), then 
ak acı, bk = 2p-1, Zk = Yk-1; Yk = ap + Arse, 


else 
dk = Yk—15 Ök = bk—1, Yk = Zk-1;Žk = bi, — A zə. 


Go to step 2. 


We can make the following remarks concerning this class of algorithms: 


1. Ay = by — ap = (—1)*81(FA — Fk-i)Ao, k > 2, where F, is the FIBONACCI 
sequence starting with Fy = Fy = 1 and Fk42 = Fku + Fk; 


2. the conditions laz, özl C [ag—1, Özil, k = 1,...,n are equivalent to certain con- 


ditions for A. Note Azı = Fk j Fiat. Then 


An < A < An, if n is odd 


C.1 
Anti < À < An, otherwise (C1) 


C.3. 


10. 
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. we get the gold section method when A = lim, sə An = (V5 — 1)/2 and then 


Aç = A*Ao. The stepnumber n is given by 
A" Ao € € 


The method is optimal from the point of view of the calculus effort; 


. the FIBONACCI search method is obtained when A = A, and then A, = 


Ao Fuk / Fu. The step number is given by 


Moreover, the method is optimal from the point of view of the error in the above 
mentioned class, because, for any A satisfying (C.1), the following inequalities 
hold: 

0 < An-2(An) < An-2(A) 


. the algorithms have linear convergence; 


. they are not stable with respect to the errors on estimating the permissible values 


of A, If A = A + 6, then 
ALÇA) = ARAN) + (—1)* 15 LAQ 


It is thus necessary to stipulate besides the error level the control of this error, 
i.e., A must be computed with an approximation less then ej (Fo); 


. for n > 25 the values of A for the gold section and for the FIBONACCI search 


method are the same in the first 10 decimal digits; 


. to counteract the shortcomings of the last two remarks, we can adopt a modified 


variant, using the FIBONACCI search method in several cycles. One cycle consists 
of running trough n — 2 steps with the standard method. Then we restart the 
process with ap = aş, bo = bn-2. After p cycles we get 


Ap(n-2)() < Ao/ F? 


. both methods assume a certain symmetry, which ensures the independence of 


the sequence of approximation intervals from the initial interval; 


the convergence conditions implies some constraints concerning A at each step: 
as upper bound at a certain step, and as lower bound at the next. 
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C.4 (o, 0)-algorithm 


We assume that we have two parameters. We expect that the conditions for the 
decrease of the approximation interval requires an upper bound for one of those para- 
meters, and a lower bound for the other, and the calculus effort remain rather identical. 

The question if we can obtain a more efficient method, in the sense to reduce the 
number of objective function evaluation. The answer is affirmative, as will be proved 
in the following. We can see in the numerical examples that the number of estimates of 
the objective function g is lowest than the FIBONACCI search method, but in addition 
a recursive calculus of an other A; is request at each step. 

We consider the following (a, 3)-algorithm, where a, 3 € İİ, 1]: 


Step 1. By a certain method, we determine an interval [ao, bo] which includes the 
solution a”. We give the error tolerance e. Set 


Ao = bo — do, Ala) = ao, A:(6) — 0 Av, 
Aə(o) = Ao — Ai(a), Az(6) = Ao — A1 (8), 
Yo = dü + Ao(a@), zo = bo — Ala), k= 0 
We compute g(yo) and g(zo). 


Step 2. If bk — a, < €, then we put «* = argmin{g(yx),¢q(z.)} and we stop the 
algorithm, else k + k +1. 


Step 3. We put Azio(a) = Akla) — Azry (a), Agge(B) = AL(9) — Anyi (6), and, if 
glyk-1) < g(zi-i), then 


Yk = ap + Agəə(a), if k is odd 
Yk ap + Aksə(0), if k is even 


? 


dk = dk-ı, Ok = Zk-1, Zh = Yk-1, { 


else 


ak = YR-1 by = by 1, Yk = 2R-1 Zk = by — Ak-ə(a), if k is odd 
a m T? | Zk = bk — Agge(B), if k is even’ 


Go to step 2. 


C.5 Properties of the (a, 9)-algorithm 


Concerning the proposed algorithm we make the following statements. 


Lemma C.1. Ifa € 9, then, for any k > 0, the inequalities 
Azizi(a) < Ares), Az(8) < Axla) 


hold. 


To prove this statement we use the first remark from the previous section. 
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Lemma C.2. For any k > 2, the equality 


holds. 


To prove this statement we use the recurrence relation for the FIBONACCI’s se- 
quence. 
The two lemmas lie to the following more general statement. 


Proposition C.1. /104] Ifa x 6, then for any k > 0, the inequalities 
min {Ag(a), An (9) } — (Pegi — 1) (8 —a)Ao < br— ar < max {Az(a), Ak(6)) 


min {Agy2(@), Appo(B)} < YR — ae, be — ə, < max {Apyo(a), Aksə(6)) 
hold. Moreover, for each strictly unimodal function q, there is a pair (ao, bo) for 


which the left equality from the first relationship holds, and another for which the 
right equality holds. 


The question is whether there are in this new class some algorithms more efficient 
than the FIBONACCI search method. First of all, that implies the stop of the new 
algorithm no later than the (n — 2)f” step, where n satisfies the inequality (C.2). 


Proposition C.2. /70// If n verifies the condition (C.2), and 


1 
=), ,a<8<\, +——— 
° 0 SOSA + ERI 


, forn even 


B=”, , An mün,“ m s?) Jorn odd 
then 
bp —ap > 0, kÜ,....mn—l 
0 < bp — dazı ) bn—2 — an—2 < An-2(An) 
hold. 


Drawing the conclusions, we can make the followings statements: 


1. the sequence of the approximation interval lenghts consists on only positive num- 
bers and it is superior bounded by a sequence of numbers which decrease to zero. 
Consequently, for n — oo, the sequence of approximations, given by this method, 
is leading to the correct solution; 


2. we have obtained a class of methods depending on a parameter, which gives the 
approximate solution bellow the error tolerance level no latter than the step, 
where the FIBONACCI search method stops; 


3. for a given function g the approximation interval lenghts depend on the position 
of the initial values ag, bọ and the ideal solution of the problem; 
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4. note that, for a = 9 = àn, we obtain the FIBONACCI search method. 


We have a whole class. The problem is to choose the optimal method. Probably, 
we reach this purpose when 5,2 — azə = 0. This equality holds when 


1 
a = àn, 9 = AQ +=, ifn is even 
Fal Fa — 1 
A A “i | if n is odd . 
B=An, a= ümu lin 18 O 


Moreover, for these values 6,1 — an-ı = 0 holds. 
The following result proves the efficiency of the proposed method. 


Theorem C.1. /96/ For a given strictly unimodal function g and an (a, 3)-algorithm 
with the parameters satisfying (C.3), and n satisfying (C.2), there is a set of (€, do, bo) 
for which the algorithm is stopped at the (n — 3)" step. 


Application. Let the problem to find the stability angle of the GEAR’s backward 
differentiation formulae. We have use the p-cyclic versions with n = 10 for both and 
(a, 3)-algorithms. The initial interval is [0,7], and the error tolerance, e = 1075. The 
results are presented in Table 


Table C.1: Approximative values of the stability angles.for GEAR’s formulae 


Formula | A-algorithm | Angle | (a, /3)-algorithm| Angle 
stepnumber | stepnumber stepnumber 
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